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ABSTHACT 

Since  widespread  use  of  the  finite  element,  method  began 
in  the  early  1960's,  much  effort  has  been  devoted  to  the 
development  of  the  method  itself,  while  only  recently  has 
there  been  any  research  directed  at  minimizing  the  discreti- 
zation error  by  a  proper  selection  of  -he  element  grid.  This 
paper  examines  some  recently  proposed  grid  optimization 
techniques  and  applies  them  to  some  one-  and  two-dimensional 
linear  self-adjoint  boundary  value  problems.  Guidelines 
requiring  minimal  software  modification  are  recommended  to 
assist  the  analyst  in  obtaining  improved  finite  element 
solutions. 
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The  critical  concern  in  any  finite  element  analysis  is 
the  adequacy  of  the  selected  mesh  to  provide  reliable  solu- 
tion results  within  seme  reascnabla  cost  range.  The  goal  of 
finite  element  grid  optimization  -hen  becomes  one  of 
obtaining  maximum  solution  accuracy  for  a  prescribed  anal- 
ysis cost.  While  this  objective  is  generally  net  realized  in 
today's  widespread  use  of  finite  eiament  analysis,  the  effi- 
cient computation  of  solutions  with  optimal  accuracy  will 
become  paramount  to  the  engineer  as  finite  element  methods 
are  applied  to  increasingly  difficult  dynamic,  nonlinear, 
and    evolution   problems. 

A.       HISTORICAL    BACKGROUND 

In  the  early  1960's,  with  the  help  of  the  high  speed 
digital  computer,  finite  element  methods  revolutionized 
problem  solving  in  engineering.  Since  that  time  the  major 
research  efforts  have  concentrated  on  expanding  the  theoret- 
ical basis  of  the  method  and  extending  its  application  in  a 
variety  of  fields.  Only  recently  has  there  been  significant 
attention  directed  at  minimizing  finite  element  solution 
errors  by  properly  defining  the  element  grid.  Early  attempts 
at  distributing  the  nodes  and  choosing  the  elements  to 
ensure  seme  degree  of  confidence  in  the  solution  accuracy 
were  predominantly  dependent  upon  :he  analyst's  engineering 
judgement  and  experience,  since  there  were  no  established 
procedures  for  accomplishing  this  objective.  Even  these 
attempts  towards  grid  optimization  have  become  largely 
ignored  with  the  advent  of  automatic  mesh  generators,  which 
have        drastically        reduced        preprocessing        costs        while 


accomplishing  little  in  improving  solution  accuracy.  These 
programs  automatically  construct  the  element  grid,  usually 
in  a  uniform  manner,  after  the  user  merely  defines  the 
problem  and  specifies  the  number  of  elements  to  be  used.  To 
establish  convergence  and  achieve  the  desired  solution  accu- 
racy, the  user  simply  repeats  the  analysis  using  a  finer 
mesh  of  uniformly  distributed  elements  while  monitoring  such 
convergence  indicators  as  successive  solution  values  at 
common  nodal  points  or  the  asymptotic  increase  in  the  energy 
content  of  the  mesh.  The  often  disastrous  flaw  in  such  a 
practice  is  that  for  nearly  degenerate  problems  which 
exhibit  very  large  gradients,  the  asymptotic  range  is  only 
entered  for  an  extremely  large  number  of  degrees-of- freedom, 
often  beyond  computer  limitations  [Ref.  1  ]•  In  this  case, 
uniform  mesh  refinement  may  produce  apparent  convergence, 
when  in  fact  the  solution  accuracy  is  poor.  Therefore,  the 
need  for  a  finite  element  grid  optimization  procedure  is 
clearly    evident. 

The  first  formal  attempts  at  finite  element  grid  optimi- 
zation did  not  begin  until  the  early  1970fs.  These  early 
approaches  involved  the  inclusion  of  the  nodal  coordinates 
as  dependent  variables  in  the  minimization  of  the  potential 
energy   functional  [Ref.    2],  Unfortunately,      the   resulting 

system  of  equations  is  highly  nonlinear  and  the  computa- 
tional effort  involved  in  its  solution  is  so  great  that 
similar  accuracy  can  be  obtained  at  a  fraction  of  the 
expense,  simply  by  employing  a  very  fine  mesh.  Clearly,  this 
method  does  not  approach  the  finite  element  grid  optimiza- 
tion goal  of  achieving  a  specified  solution  accuracy  for  a 
minimum  analysis  cost.  For  this  reason,  virtually  all  of  the 
grid  optimization  techniques  since  developed  are  based  on  a 
"near-optimum"  strategy  whereby  nearly-optimal  solution 
results  are  obtained  without  the  computational  inefficiency 
of   a    numerical    optimization      analysis.      The   growing   emphasis 
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has  been  on  adaptivity,  a  procedure  for  efficient  construc- 
tion of  near-optimum  grids  by  the  iterative  application  of 
some  criterion,  based  on  data  already  computed  from  the 
solution  for  a  previous  grid.  This  procedure  is  far  more 
efficient  than  the  conventional  approach  of  repeating  the 
analysis        using  successively        finer  uniform        grids. 

Experimental  self-adaptive  finite  element  codes  have 
recently  been  developed.  Starting  from  the  user's  initial 
idealization,  these  programs  automatically  generate  a  near- 
optimum   grid  and  solve  the    resulting    equations. 

B.       IHVESTIGATIVE   APPROACH 

In  undertaking  any  numerical  optimization  task,  the 
analyst  must  first  define  the  objective  along  with  any 
constraints  to  be  imposed  upon  the  objective  variables;  and 
finally  a  numerical  algorithm  must  be  prescribed  to  perform 
tha  optimization,  preferably  one  which  will  do  so  effi- 
ciently for  the  particular  problem.  Since  the  term  "optimum" 
most  often  refers  to  a  solution  obtained  by  mathematical 
programming,  which  is  very  inefficient  for  grid  optimiza- 
tion, a  near-optimum  grid  obtained  by  application  of  an 
adaptive  procedure,  henceforth  will  be  termed  an  optimum 
grid.  However,  before  such  a  grid  can  be  determined  to 
satisfy  the  stated  objective  of  obtaining  maximum  solution 
accuracy  for  a  prescribed  cost,  the  terms  "accuracy"  and 
"cost"  must  be  defined;  but,  more  importantly,  the  optimiza- 
tion goals  must  be  specified.  This  is  critical  because  grid 
optimization  can  oe  implemented  in  various  forms  depending 
upon  the  optimization  goals,  which  will,  in  general,  be 
determined  by  the  original  purpose  for  performing  the  finite 
element  analysis  [Ref.  3  ]•  For  example,  if  the  purpose  of 
the  analysis  is  to  evaluate  a  local  quantity,  such  as  the 
maximum      value      of    the      deoendent      variable      or   one      of      its 
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derivatives,  then  the  nodal  distribution  should  be  densest 
in  the  region  where  that,  maximum  is  determined.  If,  or.  the 
other  hand,  the  interest  is  on  an  integral  quantity  of  the 
dependent  variable  ever  a  region  of  the  domain,  then  the 
nodes  should  be  assigned  to  achieve  a  uniform  distribution 
of  the  error  over  that  region.  For  the  purpose  of  this 
investigation,  attention  will  be  focused  on  the  three  finite 
element  resultants  with  the  most  significance  in  solid 
mechanics  and  other  fields  in  which  finite  element  analysis 
is  employed:  the  maximum  value  of  the  dependent  variable,  or 
solution;  the  maximum  value  of  the  gradient  of  the  solution; 
and   the    integral  of    the    solution   over   the   domain. 

In  order  to  define  the  solution  accuracy,  it  will  be 
necessary  to  compare  the  error  in  the  solution  resultant 
obtained  using  an  optimal  grid  to  the  error  obtained  using 
some  baseline  grid  with  the  same  number  of  degrees-cf- 
freedem.  For  convenience,  the  reference  grid  chosen  will  be 
a  uniform  grid,  or  one  with  all  elements  of  the  same  order 
and  approximately  the  same  size,  with  the  understanding  that 
no  knowledgeable  analyst  would  atteapt  to  use  such  a  grid  in 
the    solution  of    a   prcblem   with    large    gradients. 

The  definition  of  analysis  cost  will  be  greatly  simpli- 
fied by  making  the  assumption  that  it  is  directly  propor- 
tional to  ths  number  of  degrees-of- freedom  used  to  obtain 
the  solution.  In  reality  the  cost  depends  on  many  factors, 
some        of        which        are  very        difficult        to        quantify. 

Understandably,  the  number  of  degrees-of-freedom  is  not  the 
sole  measure  of  computational  costs,  but  it  appears  to  be  an 
adequate  measure  of  preprocessing  and  postprocessing  costs 
which  often  account  for  the  major  portion  of  the  total 
analysis. 

This  investigation  will  concentrate  on  the  use  of  finite 
element  grid  optimization  methods  for  solving  problems  of 
structural    mechanics.  While    this      area    has      dominated  the 


literature  on  the  subject,  the  techniques  presented  herein 
extend  equally  as  well  to  any  field  for  which  variational 
principles   apply. 

There  are  two  key  questions  which  must  be  answered  prior 
to  the  adaptive  application  of  finite  element  grid 
optimization: 

(1)  What  criterion,  based  on  the  results  of  an  initial 
finite  element  analysis,  should  be  used  to  indicate 
regions   where    the   original   grid   is   inadequate    ? 

(2)  What    method   of    grid   refinement    should   be   employed   ? 


Considerable  attention  will    be  devoted  to   these   questions   in 
the    next   two  chapters. 

C.       OBJECTIVES 

The   objectives    of  this    investigation   are: 

(1)  To  examine  some  recently  developed  grid  optimization 
techniques  which  have  reached  the  state  of  practical 
applicat io  n. 

(2)  To  implement  some  of  these  techniques  in  the  solu- 
tion of  some  one-  and  two-diasnsional  linear  self- 
adjoint   boundary-value   problems. 

(3)  To  draw  a  comparison  among  t.asse  applied  techniques 
in  terms  of  solution  accuracy,  analysis  cost,  and 
ease   of  implementation. 

(4)  To  recommend  some  guidelines  to  assist  the  anaiys- 
in  obtaining  optimal  finite  element  solutions 
employing  currently  available  or  easily  amendable 
software. 
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II.    CRITERIA    FOR    GRID    REFINEMENT 

The  primary  theoretical  concern  in  the  application  of 
adaptive  grid  optimization  is  the  selection  of  the  refine- 
ment criterion.  In  ether  words,  one  must  decide  upon  which 
solution  parameter,  obtained  from  an  initial  idealization, 
may  most  appropriately  be  used  to  give  some  indication  as  to 
where  the  initial  grid  is  inadequate  and  thus  needs  refine- 
ment. There  are  several  competing  proposals  concerning  the 
most  appropriate  choice  of  a  refinement  criterion.  In 
reality,  the  decision  must  be  based  upon  such  factors  as  the 
type  of  problem  being  solved,  which  criterion  is  most  prac- 
tically implemented,  and  whether  accuracy  is  desired 
locally,  globally,  pcintwise,  or  with  respect  to  an  integral 
quantity.  The  following  are  some  of  the  more  practical 
refinement   criteria    used   in    grid    optimization   at    present. 

A.       SOLUTION  PARAMETER   7ARIATIOH 

The  most  direct,  computationally  inexpensive,  and 
earliest  proposed  indication  of  where  an  element  grid 
requires  refinement  is  a  measure  of  the  variation  of  some 
solution  parameter  over  the  domain.  It  is  only  logical  that 
a  piecewise  polynomial  approximation  would  experience  the 
most  difficulty  in  modeling  the  desired  response  in  those 
regions  where  the  solution  or  its  resultants  were  either  not 
smooth  or  were  characterized  by  large  gradients.  Therefore, 
the  basis  cf  this  criterion  is  to  refine  the  mesh  in  those 
areas  where  a  solution  parameter  varies  rapidly,  with  the 
implication  that  the  optimum  mesh  is  one  for  which  the  solu- 
tion parameter  variation  over  each  element  is  uniform 
throughout      the      domain.         The   first      consideration      in      the 
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application  of  such  a  criterion  is  to  find  a  scheme  for 
distributing  the  nodes  to  achieve  such  a  condition.  For 
one-dimensional  problems  the  task  is  trivial,  but  one  way  to 
ensure  equal  variation  over  each  element  in  higher  dimen- 
sions is  to  position  the  nodes  along  equidistant  contours  of 
the      chosen    parameter.  This      is      precisely   the      procedure 

prescribed  for  a  practical  optimization  technique  known  as 
contouring.  The  other  consideration  is  the  determination  of 
which  solution  parameter  is  to  be  usad.  In  fact,  several 
solution  parameters  have  been  found  to  work  quite  well 
[Ref.  4],  but  the  most  commonly  used  and  the  one  that  is 
consistent  with  the  potential  energy  minimization  formula" 
tion    is      the  strain    energy      density    [Ref.    2].  Because   its 

employment  requires  only  minor  software  changes  and  it  has 
been  found  to  produce  excellent  results,  this  refinement 
criterion  was  used  extensively  throughout  the  course  of  this 
research. 

B.       GRID    ITERATION 

Another,  rather  basic  but  less  direct,  method  of 
locating  regions  of  the  mesh  to  be  refined  is  known  as  grid 
iteration,  which  can  be  implemented  in  one  of  two  ways.  An 
initial  coarse  grid  analysis  may  be  repeated  using  either  a 
finer  or  a  higher  order  mesh,  and  a  comparison  of  the  resul- 
tants of  interest  between  the  two  solutions  will  identify 
those  areas  of  the  domain  where  refinement  is  most  effec- 
tive. Another  approach  is  based  on  the  assumption  that  the 
greatest  benefit  is  to  be  gained  by  refinement  in  those 
regions  where  a  small  perturbation,  like  the  introduction  of 
a  single  additional  degree- of-freedom ,  produces  the  greatest 
change  in  the  solution  or  one  of  its  resultant  parameters. 
Since  additional  degrees-of-f reedom  would  be  expected  to 
produce      the  greatest     change      in      those   regions      where      the 
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desired  response  varied  most  rapidly,  refinements  based  on 
this  method  provide  results  very  similar  to  those  obtained 
using  the  solution  parameter  variation  criterion  already 
discussed.  The  grid  iteration  method  may  at  first  appear  to 
be  more  computationally  expensive,  but  it  was  developed  to 
be  most  efficiently  implsmented  employing  a  special  family 
of  elements.  These  elements,  called  "hierarchical",  possess 
some  very  desirable  properties  for  this  application  and  will 
be    discussed  in    the    next   chapter. 

C.       EIEHENT    RESIDOALS 

Tne  major  drawback  with  refinement  criteria  based  on 
solution  parameter  variations  is  that  their  applicability 
appears  limited  to  elastostatic  problems.  For  this  reason, 
several  investigators  have  recently  developed  grid  refine- 
ment criteria  based  on  element  residuals,  which  appear  prom- 
ising for  application  to  all  types  of  finite  element 
problems,  including  nonlinear  analysis.  The  reason  for  this 
is  that  the  residual  has  essentially  the  same  meaning 
regardless  of  the  problem  type  [Ref.  5].  For  example, 
consider   the  governing  differential    equation, 

D    [    u     ]    -    f    =    0 
defined   en   seme    domain,    where    D   [     ]    is   a   linear    or   nonlinear 
differential  operator,      and   the   dependent    variable    u   and   the 
non-homogeneous    term    f  are    both   functions   of    the    independent 
variables.  Let   the     finite    element      approximation   to      the 

solution  of  the  differential  equation  be  u  -  u.  Applying  the 
differential  operator  to  the  approximation  gives  rise  to  the 
residual,    which    is    defined    as 

R    =    D    [     u    ]    -    f 
and    is  not   identically  zero      unless   the    finite   element    solu- 
tion   is   exact. 
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The  key  to  using  the  residual  as  a  criterion  for  deter- 
mining regions  of  the  domain  where  refinement  is  necessary 
is  the  local  residual  on  the  element  level,  which  indicates 
the  contribution  of  the  element  to  the  total  error  of  the 
finite  element  approximation.  Sines  the  rasidual  is  a  point- 
wise  quantity,  the  useful  measure  of  the  element  error 
contribution  is  a  norm  of  the  element  residual,  or  the  inte- 
gral over  the  element  of  the  product  of  the  residual  and 
some  weight  function.  The  integration  is  nost  readily 
performed  using  numerical  quadrature.  The  grid  optimization 
strategy  then  becomes  one  of  refining  the  mesh  so  as  to 
equi-distribute  the  element  residual  norms,  by  forcing  them 
to   become   smaller  and  more    uniform   over   the   domain. 

There  are  some  drawbacks  to  the  element  rasidual  refine- 
ment criterion  which  have  not  yet  been  fully  resolved,  such 
as  appropriate  selections  of  the  rasidual  norm  and  the 
weight  function,  and  in  the  computation  of  the  residual 
itself.  While  the  evaluation  of  the  rasidual  presents  no 
particular  difficulties  in  the  interior  of  the  element,  it 
is  rarely  determinable  at  the  edges.  The  reason  for  this  is 
that  in  formulating  the  finite  element  approximation  the 
element  shape  functions  are  defined  so  as  to  provide  only 
that  degree  of  continuity  required  to  adequately  model  the 
physical  problem;  the  most  frequent  consequence  being  that 
D  [u]  is  undefined  along  the  interelement  boundaries. 
Unfortunately,  this  singularity  cannot  be  ignored  and  a  more 
complicated  analysis  must  be  applied  in  order  to  bound  the 
residual   contributions  at   these   boundaries   [Ref.    6]. 

D.       A-POSTEBIORI    2RR0R   ESTIMATES 

A  sophisticated  extension  of  the  element  residual 
criterion  is  one  based  on  computable  error  estimates  from  an 
initial   finite    element   analysis.  This   utilizes    the   energy 
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norm  of  the  residual,  in  which  case  the  weight  function  is 
the  residual  itself.  Research  in  reliable  error  estimates 
was  pioneered  by  Babuska  [Ref.  7,  83  for  linear  quadrilat- 
eral elements  and  more  recently  by  Zienkiewicz  [Ref.  9]  for 
higher  order  elements.  The  major  difference  from  the 
residual  criterion  is  that  instead  of  equi-distributing  the 
element  residual  norms  over  the  domain,  they  are  normalized 
to  compute  error  indicators  for  the  elements,  which  are  in 
turn  used  to  compute  reliable  poiatwise  error  estimates  for 
the  solution  as  well  as  the  energy  error  over  the  domain. 
These  quantities  are  of  primary  importance  because  they 
provide  not  oaly  an  indication  of  where  additional  refine- 
ment is  most  effective,  but  also  a  measure  of  the  quality  of 
the  mesh  to  determine  whether  additional  refinement  is 
necessary  [Ref.  9].  The  optimization  strategy  is  to  obtain 
a  nearly  uniform  distribution  of  the  error  indicators 
throughout  the  domain,  which  corresponds  to  minimizing  the 
error  in  the  energy  norm.  The  refinement  procedure  may  prog- 
ress until  all  the  lccal  errors  are  within  some  prespecified 
tolerance.  While  the  practical  utility  of  such  a  refinement 
criterion  is  obvious,  the  mathematical  development  and  the 
algorithms  involved  are  rather  complicated.  However,  the 
process  is  not  computationally  expansive,  and  there  now 
exists  a  prototype  self-adaptive  finite  element  code,  FEARS, 
which   implements  this  refinement   criterion   [Ref.     10]. 
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III.     METHODS    OF    GRID    BEFIHEMENT 

Once  it  has  been  determined  where  the  initial  element. 
grid  is  inadequate  and  needs  refinement,  the  next  considera- 
tion is  how  the  idealization  in  these  areas  should  be 
improved.  The  choice  of  the  refinement  method  to  be  employed 
may  well  be  a  more  important  decision  than  the  selection  of 
one  of  the  refinement  criteria  previously  discussed,  since 
at  least  one  investigator  has  observed  that  for  a  particular 
method  of  grid  refinement,  the  various  refinement  criteria 
produce   essentially    the   same   solution    results   [Bef.    11]. 

Grid  refinement  is  the  process  of  introducing  additional 
degrees-of- freedom  into  selected  regions  of  the  finite 
element   grid,   and  may   be   performed  by    one   of    three    methods: 

(1)  The  polynomial  degree  of  the  elements  remains  fixed, 
usually  at  a  low  value,  while  the  size  of  the 
elements  is  reduced.  This  has  become  known  as  the 
h-version  since  element  size  is  commonly  denoted  by 
the   letter   h. 

(2)  The  size  of  the  elements,  usually  few  in  number, 
remains  fixed  while  the  polynomial  degree  of  the 
elements  is  increased.  This  has  become  known  as  the 
p-version  since  polynomial  order  is  commonly  denoted 
by   the  letter    p. 

(3)  The  size  of  the  elements  may  be  reduced  concurrently 
with  an  increase  in  their  polynomial  order.  This  is 
known  as  the  combined  n-  and  p-version  of  the  finite 
element   method  . 
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A.       CONVERGENCE    OF    GEID    REFINEMENT 

It  is  well  known  that  ths  finite  element  method 
converges  with  an  increasing  number  of  degr ees-of-freedcra; 
in  fact,  this  is  the  justification  for  its  development. 
Therefore,  the  appropriate  measure  of  the  effectiveness  of  a 
particular  grid  refinement  method  should  be  the  associated 
rate  of  convergence,  which  generally  will  be  affected  by  the 
smoothness  of  the  approximated  function  over  the  subdomain 
of  interest.  It  has  been  demonstrated  that  when  the  refine- 
ment is  performed  uniformly  over  the  domain,  the  associated 
rate  of  convergence  is  asymptotic,  provided  the  number  of 
degrees-of-freedom      is      sufficiently    large     [Ref.    1  ].  The 

asymptotic  rate  of  convergence  is  often  measured  as  the 
slope  of  the  error  versus  cost  curve  in  the  linear,  or 
asymptotic,  range  when  plotted  on  a  log-log  scale.  In 
performing  such  an  error  analysis  for  the  displacement 
formulation  of  the  finite  element  method,  the  error  is 
usually  the  relative  strain  energy  error,  approximated  by 
the  energy  norm,  and  the  cost  is  assumed  to  be  some  simple 
function  of  the  average  element  size  or  the  number  of 
degrees-cf- freedom      [Ref.    12:    p.    726].  Only  in      the      past 

several  years  has  there  been  any  significant  research 
comparing  the  relative  merits  of  the  different  methods  of 
grid  refinement.  Since  the  solutions  of  elliptic  boundary 
value  problems  are  usually  very  smooth  over  convex  domains 
except  in  the  vicinity  of  corners,  most  of  this  research  has 
focused  on  solutions  exhibiting  singularities,  which 
severely  hinder  the  rate  of  convergence,  as  in  problems  of 
fracture  mechanics  and  in  domains  with  re-entrant  corners 
[Ref.    I,    13,    14,    15]. 

In  order  for  a  finite  element  analysis  to  be  both  effi- 
cient and  reliable,  the  asymptotic  convergence  range  should 
he      entered     for      as    few      degrees-of-freedom      as      reasonably 
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possible.  In  general  the  p-version  satisfies  this  require- 
ment better  than  the  h-version.  While  it  has  been  estab- 
lished that  p-ccnvergence  will  necessarily  occur  whenever 
h-convergence  occurs,  the  converse  is  not  true.  For  example, 
if  the  h-version  using  a  uniform  grid  of  linear  elements  is 
applied  to  a  nearly  degenerate  problem,  the  number  of 
degrees-of- freed  cm  required  for  entry  into  the  asymptotic 
range  may  be  beyond  the  computer's  round-off  limitations,  in 
which  case  convergence  will  not  occur  unless  the  polynomial 
order  is  increased  [ Ref .  1].  Numerical  experiments  on  such 
problems  clearly  indicate  that  the  p-version  requires 
considerably  fewer  degrees-of-freeiom  than  the  h-version  to 
achieve  the  same  degree  of  accuracy.  Recent  analyses 
(Ref.  1,  13]  of  the  asymptotic  rata  of  convergence  in  energy 
for  non-smooth  solutions,  using  uniform  refinement  with 
sufficiently  high  numbers  of  degrees-of -freedom,  have  demon- 
strated that  the  p-version  cannot  have  a  slower  rate  of 
convergence  than  the  h-version.  Furthermore,  if  the  singu- 
larity is  confined  to  element  boundaries,  as  is  usually  the 
case,  the  error  for  p-method  is  inversely  proportional  to 
the  number  of  degrees-of- freedom,  whereas  the  error  is 
inversely  proportional  to  the  square  root  of  the  number  of 
degrees-of-freedom  in  the  h-version.  In  other  words,  for 
this  special  class  of  problem,  the  rate  of  convergence  for 
the  p-version  is  twice  that  for  the  h-version,  which  is  due 
primarily  to  the  ability  of  higher  order  polynomials  to 
"absorb"  singularities  occurring  at  the  element  boundaries. 
This  implies,  at  least  for  this  type  of  problem,  that  in 
order  to  minimize  the  error  for  a  specified  number  of 
degrees-of-freedom,  the  best  strategy  is  not  to  subdivide 
the  domain  uniformly,  but  to  use  instead  a  single  element  of 
increasing   polynomial  order    [Ref.    15]. 
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Since  it  is  unlikely  that  one  would  attempt  to  solve 
such  a  problem  using  uniformly  finer  grids,  a  more  useful 
comparison  between  the  convergence  rates  of  the  two  versions 
would  be  based  on  adaptive  refinement  employing  one  of  the 
solution-based  criteria  discussed  in  the  previous  chapter. 
It  so  happens  that  the  h-version,  when  used  with  optimally 
refined  meshes,  can  have  a  higher  convergence  rate  than  the 
uniformly  distributed  p-version,  provided  that  the  element 
order  is  sufficiently  high.  However,  the  p-version  can  also 
be  employed  with  an  optimal  refinement  criterion.  While 
there  are  yet  no  proven  theorems  concerning  the  convergence 
rates  for  non-uniform  refinement,  obtaining  the  desired 
solution  accuracy  with  optimal  p- distributions  appears  to  be 
much  less  sensitive  to  the  particular  choice  of  the  elements 
to    be   refined  than    with   optimal   h-refinement   [Ref.    13]. 

It  would  seem  plausible  that  an  even  better  optimization 
strategy  would  involve  a  proper  combination  of  both  the  h- 
and  p-versicns.  It  has  been  demonstrated  for  problems  with 
corner        singularities,  that        a        proper        sequence        of 

h-ref inements  combined  conourrently  with  the  proper  sequence 
of  p-distributions  results  in  extremely  high  convergence 
rates,  conjectured  tc  be  exponential  [Ref-  15].  However, 
this  proper  combination  is  difficult  to  determine,  and  adap- 
tive refinement  based  on  the  combined  h-  and  p-versions 
poses  some  difficult  data  management  problems.  To  avoid  this 
problem  a  more  promising  approach,  proposed  by  Babuska  and 
Szabo  [Ref.  1],  employs  a  graded  mesh  in  which  the  element 
sizes  are  first  reduced  according  to  a  geometric  progression 
towards  the  singularity,  followed  by  determining  the  optimal 
p-distribution  for  those  elements  using  an  adaptive 
criterion.  However,  obtaining  the  optimal  combination  when 
employing  this  scheme  can  be  a  delicate  matter  and,  astcund- 
ingly,  the  highest  accuracy  is  achieved  when  the  polynomial 
order  of  the  elements  actually  increases  with  distance  from 
the   singularity. 
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There  are  seme  additional  advantages  of  the  p-version 
worth      mentioaing.  Because   the      p-version      employs      fewer 

elements,  there  are  lesser  preprocessing  and  postprocessing 
costs  than  foe  the  h-version.  Furthermore,  when  bandwidth 
minimization  and  sparse  matrix  solution  techniques  are  used, 
the  solution  time  for  the  p-version  is  approximately  the 
same  as  for  the  h-version  for  a  specified  number  of 
degrees-of- freedom,  and  the  p-version  appears  less  suscep- 
tible to  round-off  errors.  Finally,  the  p-version  is  simpler 
to  implement  adaptively  than  the  h-version  when  hierarchical 
elements   are  employed  [Ref.     13]. 

B.       HIERARCHICAL    FINITE    ELEMENTS 

The  hierarchical  concept  was  first  introduced  as  a 
simple  method  for  implementing  the  p-version  and  as  a 
convenient  device  fcr  imposing  boundary  continuity  between 
elements  of  different  polynomial  order  [3ef.  9].  Since  then 
a  useful  family  of  elements  based  on  the  hierarchical 
concept  has  been  developed  and  incorporated  into  COMET-X,  an 
experimental  finite  element  cede,  developed  by  Szabo,  which 
self -adaptively  employs  both  the  h-  and  p-versiens  of  the 
finite   element    method  [Ref.     14]. 

For  a  brief  description  of  the  hierarchical  concept 
consider  the  conventional  finite  element  formulation  which 
produces   the  following  system   of   equations: 


K,   xu<n)=      f(n)  (Eqn.    3.1) 


where      n      is     the  number   of    degrees-of- freedom, 
n  x    n      global   stiffness      matrix,      u^  is   the      finite   element 
approximation     of      the     exact      solution,  ar.d      f'n'    is      the 

n-component     global      lead      vector.      When        m        higher      order 
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degrees-of-freedcm  are  added  to  the  original  system  using 
conventional  refinement  methods  the  system  of  equations 
becomes: 


$n-Hn)  « 


^ 


(Eqr.    3.2) 


where   the      contributions   to   £(n+m\    and   f 


(n+m) 


from      the   refined 


elements  result  in  a  completely  different  set  of  coeffi- 
cients. If,  on  the  other  hand,  this  refinement  had  been  made 
hierarchically,    the    equations   would    become: 
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(Eqn.    3.3) 


where    K, 


and  f,„\are   the      stiffness    matrix    and      rorce   vector 


tfn)  *iiU  «n) 
from    the    original     system  of      equations    for        n      degrees-of- 

freedom  appearing  in  Equation  3.1.  However,  u'm'is  not  the 
nodal  values  of  the  finite  element  solution  for  the  m 
additional  degrees-of-fr  eedom,  but  instead  represents  the 
difference  between  those  values  and  the  pointwise  values 
obtained  from  the  lower  order  polynomial  interpolation  for 
the    unrefined  mesh    of  n   deg rees-of-f reedom. 

The  primary  advantage  of  hierarchical  elements  is  imme- 
diately observable  from  Equation  3.3.  Because  the  shape 
functions  of  an  element  of  order  p  constitute  a  subset  of 
the  shape  functions  of  an  element  of  order  p+1,  the  local 
stiffness  matrix  and  force  vector  for  each  hierarchical 
element  is  embedded  in  the  stiffness  matrices  and  force 
vectors  of  all  hierarchical  elements  of  higher  order. 
Therefore,  the  global  stiffness  matrix  K/  vand  force  7ector 
f        of      the     original  system      are      preserved,        thus      saving 
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considerable  time  and  effort  expended  on  computing  the  coef- 
ficients for  successive  refinements  [Ref.  14].  Another 
advantage  is  that  the  hierarchical  form  of  the  global  stiff- 
ness matrix  is  more  diagonally  dominant  than  the  one 
resulting  from  a  conventional  refinement,  resulting  in 
improved  conditioning  and  faster  convergence  whan  iterative 
solvers  are  employed  [Ref-  9].  Another  benefit  of  hierar- 
chical elements,  which  arises  from  the  "add-on"  nature  of 
the  nodal  variables  cf  the  higher  order  degrees-of- freedom, 
is  that  the  problem  of  maintaining  boundary  continuity 
between  elements  of  different  polynomial  order  becomes 
trivial.  Instead  of  introducing  global  constraint  equations 
for  the  higher  order  degrees-cf-freedom,  the  nodal  variables 
are  simply  set  equal  to  zero  and  condensed  out,  as  if  they 
were    zero-valued   Dirichlet    boundary    conditions  [Ref-    2]- 

There  are  two  major  drawbacks  with  hierarchical  elements 
that  have  hampered  their  widespread  acceptability.  The 
first,  which  has  already  been  mentioned,  is  that  the  nodal 
variables  for  the  higher  order  degrees-of-fr eedom  represent 
difference  values  rather  than  the  aore  easily  identifiable 
values  of  the  dependent  variable  itself.  Secondly,  when 
implementing  the  h-version  cf  the  finite  element  method, 
special  integration  rules  must  be  introduced  when  the  subdi- 
vided element  is  in  hierarchical  form  [Ref.  9]-  Of  course, 
the  latter  problem  can  be  evaded  by  using  the  p-version,  for 
which  the  hierarchical  concept  was  developed.  In  spite  of 
the  disadvantages  of  hierarchical  elements,  their  consider- 
able computational  efficiency  and  utility  for  grid  optimiza- 
tion will  certainly  result  in  their  widespread  utilization 
in    future  adaptive    finite   element   oodes. 
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17.    GRID   OPTIMIZATION    TECHNIQUES 

Once  the  analyst  has  identified  where  the  initial  grid 
needs  enrichment  and  decided  which  refinement  method  to 
employ,  he  mast  then  determine  a  systematic  procedure,  or 
algorithm,  to  perform  the  refinement  according  to  the 
criterion  selected.  The  ultimate  goal  of  such  a  procedure  is 
to  design  an  element  grid  which  meets  the  optimization 
objective  of  obtaining  maximum  solution  accuracy  for  a  spec- 
ified analysis  cost.  While  the  analyst  may  or  may  not  have 
an  indication  of  the  accuracy  of  the  solution,  he  should 
have  a  preconceived  notion  of  cost,  or  how  much  effort  he  is 
willing  to  expend  to  arrive  at  a  better  solution.  Therefore, 
with  some  knowledge  of  the  grid  optimization  techniques 
available  and  an  understanding  of  the  advantages  and  disad- 
vantages of  each,  the  analyst  can  realize  the  grid  optimiza- 
tion   goal. 

There  are  essentially  two  adaptive  grid  optimization 
strategies: 

(1)  Grid  refinement,  in  which  the  initial  analysis  is 
performed  on  a  relatively  coarse  grid,  and  new 
degrees-of -freedom  are  added  to  the  same  grid  by  the 
iterative  application  of  the  solution-based 
criterion. 

(2)  Grid  modif ica ticn,  in  which  the  initial  analysis  is 
performed  using  a  prespecified  number  of 
degrees-of-freedom,  and  the  solution-based  criterion 
is  employed  to  shift  degrees-of-freedom  from  certain 
regions  to  others.  This  may  involve  complete  grid 
redefinition  in  an  effort  to  obtain  a  near-optimum 
grid   in  a    single  cycle. 
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Much  of  the  interest  lately  has  been  in  the  development 
cf  complicated  self-adaptive  software  packages  which  mini- 
mize the  impact  of  the  user's  skill  on  the  final  solution. 
Ideally,  the  analyst  would  merely  define  the  problem  and  tha 
program  would  automatically  generata  and  analyze  the  optimum 
grid  employing  one  or  more  of  these  techniques,  possibly  at 
the    user's   cption. 

A.       MATHEMATICAL    PROGRAMMING 

No  discussion  of  grid  optimization  techniques  would  be 
complete  without  a  brief  description  of  mathematical 
programming,  not  only  because  it  is  how  grid  optimization 
was  earliest  attempted,  but  more  importantly,  it  is 
precisely  what  the  engineer  envisions  when  he  hears  the  term 
"optimization".  It  is  not  a  grid  optimization  technique,  per 
se,  but  rather  a  numerical  process  of  achieving  any  optimi- 
zation objective,  once  it  is  explicitly  defined  in  mathemat- 
ical terms.  In  solid  mechanics  the  finite  element  method  is 
a  numerical  method  for  minimizing  the  potential  energy  func- 
tional,   which  in   discretized   form   may   be    written: 


7T  *  k    uT    K   u    -    uT   f  (Eqn.    4.1) 


where:  u  is  the  global  displacement  vector 


a. 


K    is  the  global   stiffness  matrix,    and 
f   is  the   global   load   vector 


In   the   classical   finite      element  formulation,      the   potential 

energy   is    minimized    with   respect  to    the    nodal    displacements, 

which      implies      satisfaction      of  the      following      stationary 
conditions: 
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3u4     =0  (i   =    1,2,.. .,n)  (Eqn.    4.2) 


i 


where  n    is   the    number     of  degrees-of- freedom.      This   leads   to 
the    very    familiar  system   of    linear  aquations: 


{   U   -    f    ■   0  (Eqn.     4.3) 
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However,  since  K  and  f  are  functions  of  the  nodal  coordi- 
nates, then  the  potential  energy  could  be  minimized  with 
respect  to  the  nodal  coordinates  as  well.  This  would  require 
satisfaction  of  the  following  additional  stationary 
conditions: 


3TT 

3xj      =0  (j    =     1,2,--.,  m)  (Eqn.    4.4) 


where  m  is  the  number  of  nodal  coordinates,  x..  This  differ- 
entiation leads  to  the  less  familiar  system  of  non-linear 
equations: 

t  3&  3fT 

k  u'  -e*-!!  -  -^-  u  =  0  (J   =    1,2,.. .,m)  (Eqn.    4.5) 

"*  %     3x.   /\J       3x-   a, 

This,  then,  is  the  mathematical  statement  of  the  grid  opti- 
mization problem  for  the  elastostatic  case.  The  nodal 
displacement  variables  may  be  eliminated  by  minimizing  the 
potential  energy  with  respect  to  the  nodal  coordinates  only, 
subject  to  the  implicit  constraint  that  Equation  4.3  is 
always  satisfied  [Ref.  4].  Unfortunately  this  does  not  help 
much  because  the  objective  function  is  still  nonlinear, 
rendering  most  numerical  optimization  algorithms  inefficient 
and    unreliable.    The    difficulty    is   even    further  compounded   by 
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the  requirement  that  the  nodal  variables  be  subject  to  side 
constraints  in  order  to  maintain  the  defined  boundary  of  the 
domain  and  to  ensure  that  the  elements  neither  distort 
excessively  nor  overlap  one  another.  For  all  except  the 
simplest  of  problems,  these  constraints  may  be  even  more 
severely  nonlinear  than  the  objective  function,  resulting  in 
the  analysis  becoming  prohibitively  expensive  [Ref-  2].  For 
this  reason,  mathematical  programming  in  finite  element  grid 
optimization  has  been  abandoned  in  favor  of  some  equally 
reliable,  yet  far  mere  computationally  efficient  grid  opti- 
mization techniques.  However,  these  early  efforts  with 
mathematical  programming  were  not  totally  in  vain  because 
they    gave   rise    to   the  contouring   techniques. 

B.       COHTOOHING 

Since  mathematical  programming  is  infeasible  for  grid 
optimization,  further  investigations  were  conducted  to 
suggest  some  guidelines  to  enable  the  analyst  to  construct  a 
grid  with  similar  topological  features  of  the  numerically 
optimized  grid  without  the  computational  effort  involved. 
Turcke  [Ref.  4],  in  employing  mathematical  programming  in 
the  solution  of  some  simple  two-dimensional  eiastcstatic 
problems,  observed  that  there  was  a  very  definite  element 
pattern  common  among  problems  involving  high  strain  gradi- 
ents and  that  the  nodes  of  the  numerically  optimized  grid 
generally  tended  to  be  aligned  along  contours  of  some 
response  function  being  modeled.  Consequently,  in  performing 
analyses  on  grids  whose  construction  was  based  on  contours 
derived  from  an  initial  analysis,  it  was  determined  that  the 
following  provided  grid  characteristics  in  regions  of  high 
strain  gradients  similar  to  the  numerically  optimized  grid, 
but    at    a   fraction   of  the   computational   expense: 
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•  contours  of    displacement 

•  contours  of    maximum   principal   stress 

•  contours  of    maximum   shear   stress 

•  contours  of    strain  energy   density    (isoenergetics) 

•  principal   stress    trajectories    (isostatics) 

Since  the  strain  energy  density  is  the  response  which  is 
consistent  with  the  principle  of  minimum  potential  energy, 
isoenergetics  are  the  most  commonly  used  contours  along 
which  element  edges  are  aligned  [Ref.  4  ]•  However,  there 
still  remains  the  question  of  how  to  position  the  nodes 
along  the  contours.  For  this  reason,  isostatics  have  beccme 
increasingly  popular  because  the  principal  stress  trajecto- 
ries form  a  "flow  net"  of  orthogonal  curves  which  can  guide 
the   analyst    in    the    layout   of   the   elements   [Ref.     16]. 

Since  contouring  involves  the  redefinition  of  the  grid, 
as  opposed  to  a  grid  refinement,  its  primary  advantage  is 
that  the  enriched  mesh  is  not  constrained  to  the  element 
configuration  of  the  previous  mesh.  Therefore,  there  is  no 
limit  to  the  amount  of  enrichment  per  cycle  which  can  be 
performed  and  it  is  conceivable  that  an  optimum  mesh  could 
be  generated  in  a  single  cycle  [Ref.  2].  However,  while  the 
computational  cost  of  repeated  analyses  is  reduced,  the 
preprocessing  costs  involved  in  constructing  the  contours 
and  redefining  the  mesh  can  become  quite  high,  especially  if 
the  contours  are  complex.  Algorithms  for  performing  these 
tasks  in  two-dimensional  domains  have  been  proposed 
[Ref.  Ur  11  ],  but  they  are  not  extendable  to  three- 
dimensional  problems.  The  major  obstacle  for  two-  and 
three-dimensional  domains  is  that  it  is  often  difficult  to 
constrain  the  element  edges  to  the  contours  without  the 
elements  becoming  elongated  or  distorted  to  the  degree  that 
numerical  inaccuracies  result.  Another  difficulty,  not 
addressed   in  the   literature,      is      how   the   contour    increments 
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should  be  selected  when  the  response  function  is  non- 
monotonic over    the    domain. 

C.       SELECTIVE  REFINEMENT 

The  most  commonly  employed  grid  optimization  technique 
is  that  of  selective  refinement.  As  its  name  implies, 
selected  elements  from  a  given  mesh  are  enriched  while  the 
original  element  grid  remains  sssentially  intact.  The 
elements  selected  for  refinement  are  determined  by  the  iter- 
ative application  of  the  solution-based  criterion  to  indi- 
cate which  elements  contribute  most  to  the  solution  error. 
The  refinement  can  be  performed  by  either  the  h-versicn  or 
the  p-version,  or  even  the  combined  version  if  so  desired, 
but  the  choice  is  most  often  predetermined  by  the  capabili- 
ties of  the  available  preprocessor.  Since  the  addition  of 
new  degrees-of -freedom  over  several  iterations  can  quickly 
enlarge  the  problem,  it  is  advisable  to  perform  the  initial 
analysis  with  a  reasonably  coarse  grid  of  optimally  shaped 
elements,  that  is  nearly  square  quadrilaterals  or  nearly 
equilateral  triangles.  This  is  especially  important  in  the 
h-version  where  it  is  desirable  to  prevent  the  successive 
subdivision  of  elements  from  producing  elongated  new 
elements.  One  refinement  technique  which  will  ensure  this  is 
the  so  called  "father-to-four  sens"  subdivision  scheme  in 
which  a  single  quadrilateral  or  triangular  element  is 
replaced  by  four  new  ones  by  adding  and  connecting  midside 
nodes  on  the  edges  of  the  original  element  as  shewn  in 
Figure   4.1.  The    major   difficulty    in      selective   refinement 

arises  when  the  addition  of  a  node  along  an  edge  of  the 
element  to  be  subdivided  creates  a  higher  polynomial  ordered 
edge  for  an  adjacent  element  which  is  not  tc  be  subdivided. 
There  results  an  incompatibility  in  the  interpolation  of  the 
dependent   variable    along    this   interelement   boundary.    Such   is 
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Figure  4.1         Seme   h-Version   Subdivision   Schemes. 

the  case  in  the  h-version  scheme  of  Figure  4.1  and  it  also 
arises  in  the  p- version  when  two  elsments  of  different  poiy- 
ncmial  order  share  a  common  edgs.  When  -his  situation 
occurs,  the  additional  de grees-of- freedom  do  not  actually 
represent  degrees-of-freed  om  at  all  because  they  must  be 
numerically  constrained  to  the  polynomial  interpolant  of  the 
lower  order.  Such  constraints  are  usually  imposed  in  one  of 
three  ways:  global  constraint  aqua-ions  may  be  written;  the 
constraints  may  be  incorporated  in  the  elemental  basis;  or 
hierarchical  forms  may  be  used  with  the  excess  degrees-of- 
freedem   simply    set    to   zero    and   condensed   out   [Ref.    2]. 

There  are  some  ether  selective  refinement  techniques 
which  do  not  require  any  major  software  modifications.  In 
the  h-version,  the  continuity  problem  may  be  circumvented  by 
employing  any  of  the  coarse-to-fine  mesh  transition  schemes 
for  which  all  of  the  element  edges  remain  of  the  same  poly- 
nomial order  [Ref.  17:  p.  210].  HDwe'/er,  it  is  impossible 
to      employ    these      schemes      without      permitting   some      element 
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distortion,  and  the  refinement  must  nearly  always  be 
performed  interactively  rather  than  automatically.  For  the 
p-version,  interelement  continuity  can  be  easily  ensured  by 
employing  vac  i  able-  node  d  isoparametric  elements,  which 
permit  a  single  element  to  possess  edges  of  differing  poly- 
nomial  orders  [Ref.     17:    p.     125]. 

The  analyst  must  also  exercise  care  when  adding  new 
nodes  to  the  boundary  of  the  domain  to  ensure  that  the 
appropriate  boundary  conditions  are  determined  and  applied. 
Furthermore,  if  the  boundary  is  curved,  the  coordinates  of 
the  new  node  should  be  computed  such  that  it  is  placed  on 
the  actual  boundary  and  not  necessarily  on  the  edge  of  the 
element    being  refined  [Ref.    2]. 

The  important  advantage  of  the  selective  refinement 
technique  is  that  once  an  appropriate  refinement  criterion 
has  been  determined,  selecting  candidate  elements  for 
refinement      in       each    cycle      becomes      straightforward.  The 

refinement  can  then  be  continued  indefinitely  to  achieve 
very  high  accuracy,  but  because  the  solution  phase  is 
repeated  for  each  cycle,  it  is  desirable  to  hold  the  number 
of  cycles  to  a  minimum.  Because  the  nodes  from  the  previous 
mesh  remain  fixed  for  each  cycle,  selective  refinement  is 
ideally  suited  for  iterative  solution  methods.  The  solution 
values  obtained  from  the  previous  cycle,  combined  with 
interpolated  values  for  the  new  degrees-of -freedom,  provide 
an   excellent   initial    guess    for   the   next    cycle   [Ref.    2]. 

The  major  disadvantage  is  that  the  limited  amount  of 
refinement  which  can  be  performed  in  each  cycle  may  necessi- 
tate several  cycles  tc  obtain  an  optimum  grid.  In  addition, 
if  new  degrees-of-f reedom  require  interelement  continuity 
constraints,  data  management  can  become  cumbersome  unless 
the    constraint    is   performed    hierarchically. 
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D.       SOBDOMAIN   ISOLATION 

One  of  the  obvious  disadvantages  of  the  selective 
refinement  technique  is  that  the  solution  must  be  completely 
repeated  for  each  cycle  when,  in  fact,  -he  number  of 
degr ees-cf- frsedom  added  in  each  cycle  may  be  few  in  compar- 
ison to  the  total  for  the  problem.  In  addition,  the  number 
of  elements  requiring  refinement  in  each  cycle  may  only 
account  for  a  small  portion  of  the  domain.  Although  the 
refinement  criterion  has  indicated  where  the  grid  is  inade- 
quate and  the  approximation  is  likaly  to  be  poor,  the  solu- 
tion is  repeated  in  each  cycle  for  those  nodes  where  the 
error  is  presumably  small.  Besides  the  apparent  computa- 
tional inefficiencies,  this  shortcoming  severely  restricts 
the  amount  of  refinement  which  can  be  performed  in  the 
subregicns  of  interest  since  it  is  desirable  to  confine  the 
size  of  the  problem  within  reasonable  limits.  An  alternative 
approach  is  to  reformulate  the  problem  for  those  subregions 
where  refinement  is  necessary  and  to  accept  the  results  of 
the  initial  analysis  as  an  adequate  solution  for  the 
remainder  of  the  domain.  The  elements  requiring  further 
refinement,  which  constitute  isolated  subdomains  of  the 
original  prcbleiD,  can  generally  be  subjected  to  signifi- 
cantly greater  refinement  than  would  otherwise  be  practical. 
The  solution  obtained  from  the  initial  analysis  is  then  used 
in  imposing  boundary  values  on  those  degrees-of-f reedcra 
located  on  the  boundaries  of  an  individual  subdcmain.  These 
can,  in  turn,  be  used  to  generate  the  boundary  conditions 
for  any  additional  boundary  degr ees-of-f reedom  introduced  by 
the   refinement    using   an  appropriate   interpolation    scheme. 

This  grid  optimization  technique,  which  the  author  terms 
"subdcmain  isolation",  has  some  further  advantages  over 
selective  refinement.  The  subdomain  aay  be  selected  arbi- 
trarily   small      such    that    excellent      results   may      be   obtained 
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with  a  single  cycle  using  uniform  refinement.  Therefore,  the 
difficulties  involving  coarse-to-fine  transition  schemes, 
element  elongation  and  interelement  continuity  can  be 
avoided.  Furthermore,  one  can  choose  as  many  subregions  for 
refinement  as  desired  without  creating  an  excessively  large 
problem. 

The  obstacle  which  may  prevent  this  technique  from  being 
readily  accepted  is  the  notion  that,  by  imposing  erroneous 
boundary  conditions  on  the  subdomains,  the  convergence  of 
the  finite  element  method  to  the  exact  solution  in  these 
regions  has  somehow  been  tampered  with.  This  aversion  may  be 
somewhat  abated  by  considering  a  simple  extension  of 
Saint-Venant *3       Principle  [  Ref .    18:    p.    33],  Although      the 

conditions  are  not  rigorously  satisfied  at  the  boundary, 
which  may  result  in  significant  changes  in  the  response 
locally,  the  affect  at  some  sufficient  distance  away  will  be 
negligible.  The  numerical  evidence  supports  this  premise. 
While  errors  in  the  boundary  values  may  somewhat  restrict 
the  accuracy  af  the  dependent  variable,  great  improvements 
can  be  realized  in  the  accuracy  of  its  gradients,  which  is 
more  often  the  goal  of  the  optimization.  Since  the  initial 
analysis  provides  the  boundary  values  for  the  subdomains,  it 
is  desirable  that  its  solution  be  as  accurate  as  reasonably 
possible.  Fortunately,  since  subsequent  refinements  are  not 
performed  on  the  original  grid,  the  initial  analysis  may 
involve  significantly  more  degrees-cf-freedom  than  in  the 
case    of   selective  refinement. 

E.        HESH    GRADING 

The  final  grid  optimization  technique  to  be  discussed 
employs  a  mesh  for  which  the  element  sizes  are  successively 
reduced,  according  to  some  geometric  sequence,  towards  a 
selected   region      of    the  domain.        One    might   argue      that    mesh 
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grading  is  not  really  an  optimization  technique  since  it  is 
most  often  applied  on  an  "a-priori"  basis  rather  than  adap- 
tively,  and  that  it  does  not  lend  itself  well  to  the  itera- 
tive application  of  a  solution-based  criterion.  However,  the 
technique  is  simple  to  use  and  its  implementation  requires 
few  software  modifications.  Furthermore,  a  solution-based 
refinement  criterion  can  be  used  to  give  a  measure  of  the 
quality  of  the  mesh  to  indicate  whether  a  more  pronounced 
grading  may  prove  beneficial.  Depending  on  the  solution 
parameter  of  interest,  mesh  grading  can  provide  excellent 
accuracy  at  a  low  analysis  cost.  This  refinement  method  must 
therefore  be  considered  among  the  grid  optimization 
techniques. 

For  the  less  elaborate  finite  element  preprocessors, 
mesh  grading  is  often  the  only  refinement  means  available 
without  resorting  to  a  uniformly  finer  mesh  involving  many 
more  degrees-of-freedcm.  The  most  common  method  of  implemen- 
tation in  two-dimensions  is  to  first  define  the  problem 
domain  in  terms  of  a  curvilinear  quadrilateral  by  selecting 
four  keynodes  along  the  problem  boundary.  Then  the  boundary 
nodes  are  spaced  according  to  some  geometric  sequence  based 
on  the  user-provided  bias  parameters  which  determine  the 
density  of  the  nodes  towards  selected  points  on  the  four 
quadrilateral  edges.  Finally,  curves  are  generated  to 
connect  the  boundary  nodes  on  opposite  edges  of  the  quadri- 
lateral, thus  producing  a  graded  mesh.  This  process,  which 
can  also  be  extended  to  three-dimensions,  is  the  mesh  gener- 
ation scheme  employed  in  the  finite  element  code  GIFTS 
[Ref.    19]. 

The  major  disadvantage  of  mesh  grading  is  that  in  order 
to  achieve  sufficiently  small  elements  in  the  region  of 
interest,  the  elements  must  grow  successively  larger  away 
from  that  region.  This  may  be  very  undesiraole,  especially 
if   refinement  is      called    for   in   more    than   one      region   of   the 
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domain,  in  which  case  the  mesh  must  be  generated  and  graded 
by  subdcmainsr  thereby  complicating  the  data  management 
involved. 

Another   disadvantage   is    that    unless   the   domain    possesses 
some        special        geometric      symmetry,  excessive        element 


Figure   4.2        Graded  aesh    for   a   Perforated   Square   Plate. 

elongation  will  usually  result  if  a  highly  pronounced 
grading  is  rsquired.  Some  configurations  are  particularly 
well  suited  for  refinement  for  mssh  grading  such  as  the 
classical  perforated  square  plate  problem  in  solid  mechanics 
shown    in   Figure    4.2. 


37 


V-    APPLICATION    TO    OJE^  DIMENSIONAL    PROBLEMS 

Now  that  the  necessary  tools  for  performing  grid  optimi- 
zation have  been  intrcducai ,  it  is  time  to  employ  them  in  an 
attempt  to  obtain  optimal  solutions  to  some  practical  prob- 
lems in  engineering.  An  obvious  starting  point  for  such  an 
investigation  is  the  one-dimensional  boundary  value  problem. 
While  most  of  the  fruitful  research  in  grid  optimization  has 
concentrated  on  problems  of  higher  dimensions,  the  cne- 
dimensional  problem  is  a  very  convenient  device  for  studying 
finite  element  grid  optimization.  Foremost,  one-dimensional 
finite  element  models  possess  a  unique  connectivity  in  that 
adjacent  elements  meet  at  their  end  nodal  points.  Therefore, 
refinement  by  the  h-  or  p-versionsf  or  by  relocating  nodal 
points  becomes  a  trivial  task,  which  does  not  involve  any  of 
the  difficulties  so  frequently  encountered  with  higher 
dimensional  problems,  such  as  preserving  interelement  conti- 
nuity and  maintaining  optimal  element  shapes.  Furthermore, 
one- dimensional  studies  can  often  provide  valuable  insight 
to  the  solution  of  more  difficult  higher  dimensional 
problems. 

The  primary  concerns  in  the  selection  of  the  problems  to 
fce   studied   were    as    fellows: 

(1)  there  should  exist  an  analytical  solution  to  provide 
a   means  of   reliable   error   analysis; 

(2)  the  solution  and  its  resultants  should  exhibit 
sufficiently  high  gradients  so  that  the  effective- 
ness  of  the  grid  optimization    is   readily   observed. 
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Because  of  the  complexity  and  a  certain  degree  of  arbi- 
trariness involved  in  the  computation  of  element  residuals 
and  a-posteriori  error  estimates,  the  solution  parameter 
variation  is  the  refinement  criterion  of  choice.  There  are 
several  solution  parameters  which  are  easily  computed, 
requiring  minimal  software  changes  to  an  existing  finite 
element   code. 

Furthermore,  for  the  one-dimensional  investigation,  it 
was  decided  to  simplify  the  analysis  by  exploiting  the 
linear  elements.  While  it  is  granted  that  improved  solution 
accuracy  may  generally  be  obtained  by  employing  higher  order 
elements,  it  will  be  assumed  that  conclusions  based  on  the 
use  of  linear  elements  can  be  applied  as  well  to  elements  of 
higher   pclynomial   order. 

A.       ELASTIC    CABLE   PROBLEM 

Consider  an  elastic  cable  under  tension  r,  stretched 
between  two  points  a  distance  2L  apart.  If  the  cable  is 
supported  by  a  Winkler,  or  elastic,  foundation  of  modulus  k, 
and  a  concentrated  load  P  is  applied  at  the  midspan,  the 
resulting      deflection   v(x),  (0    <   x    <   L)  ,      is      as    shown     in 

Figure  5.1.  The  analytical  solution  and  finite  element 
approach    for  this   problem   are   presented    in    Appendix    A. 

The  initial  finite  element  analysis  was  performed  using 
ten  linear  elements  of  uniform  length.  From  this  initial 
analysis  the  approximate  distribution  of  the  following  solu- 
tion   resultants    was    obtained   over   the    domain: 

•  the    displacement,   v    (the   solution) 

•  the   slope,    V 

•  the    strain    energy,   U 

•  the   strain    energy  density,    SED    (dU/dx) 
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Figure  5.1        Tension   Cable    Deflection   on   Elastic   Foundation. 

Subsequent  analyses  were  performed  for  finite  element,  models 
using  the  same  number  of  elements,  but.  with  the  nodes  redis- 
tributed to  achieve  approximately  uniform  variation  of  the 
above  parameters  over  each  element.  Note  that  the  strain 
energy  refinement  criterion  produces  elements  of  identical 
strain  energy  content.  In  addition,  the  problem  was  solved 
employing  graded  element  models  of  various  adjacent  element 
length  ratios.  The  resulting  element  models  based  on  these 
refinement  criteria  are  shown  in  Figure  5.2  (a-f).  The 
graded  model  (b)  for  an  element  length  ratio  of  1.2  is 
presented  for  comparison  because  it  produced  good  overall 
solution   results. 

As  previously  mentioned,  the  solution  resultants  of 
primary  interest  are  the  maximum  displacement,  the  maximum 
slope,  and  the  integral  of  the  displacement  over  the  domain, 
because  they  represent  important  analogous  solution  results 
in  nearly  all  fields  in  which  finite  element  analysis  is 
often   performed.         The  accuracy   obtained   in   these    values    for 
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Figure   5.2        Tension   Cable    Refinements. 
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each  of  the  refined  models  is  presented  in  Table  I.  As  car- 
be  seen  in  Figure  5.2,  the  strain  energy  and  strain  energy 
density  criteria  produced  extreme  variations  of  element 
length   while  the      criteria    of   displacement    and      slope   result 


TABLE  I 

Tension  Cable 

Problem 

Solution  Results 

Problem  Paramets 

>rs :    L 
k 
T 
P 

= 

100  in 
1  psi 
1000  lb 
1000  lb 

Variation 

Refinement 

Criterion 

Percent 
v  (max) 

age  Relative 
v'  (max) 

Error 
lvdx 

Uniform 

-0.40 

0.36 

0.  12 

Graded  (1.  2) 

-0.  19 

0.  07 

0.  17 

V 

-0.  18 

0.06 

0.39 

V' 

-0.23 

0.05 

0.87 

0 

-1.03 

0.  05 

3.53 

SED 

-1.29 

0.05 

4.17 

U    (mod) 

-0.53 

0.04 

1.63 

SED  (mod) 

-0.51 

0.03 

1.48 

in  more  moderate  variations.  It  can  be  observed  in  Table  I 
that  the  more  pronounced  refinements  based  on  energy  distri- 
bution result  in  greater  accuracy  for  the  maximum  slope  but 
with  the  accompanying  severe  penalty  of  significantly  poorer 
estimations  of  the  maximum  displacement  and  the  integral 
quantity.  For  this  particular  problem  the  uniform  grid 
provides     optimal        accuracy      of        the      integral         quantity, 
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therefore  refinement  cannot  reduce  its  error.  Yet  great 
improvement  in  the  accuracy  of  the  maximum  slope  and  modest 
improvement  in  the  accuracy  of  the  maximum  displacement  can 
be  achieved  with  moderate  refinements  based  on  the  displace- 
ment   and   slope    distributions. 

One  might  assume,  and  correctly  so,  that  the  ability  of 
the  energy  refinements  to  produce  the  best  accuracy  for  the 
maximum  slope  is  due  to  the  extremsly  small  elements  which 
result  in  the  area  where  that  quantity  occurs.  Furthermore, 
it  would  be  correct  to  propose  that  the  reason  for  these 
refinements  producing  poorer  estimates  than  the  uniform 
model  for  the  other  two  quantities  of  interest  is  that  the 
excessively  large  elements  in  the  regions  of  low  gradients 
severely  overstiffen  the  model  there.  It  would  then  seem 
plausible  tc  improve  the  accuracy  for  the  maximum  displace- 
ment and  the  integral  quantity  by  redistributing  the  ncdes 
in  these  regions  to  prevent  such  excessively  large  elements. 
This  was  dene  by  arbitrarily  employing  a  grading  scheme  to 
the  three  largest  elements  to  produce  the  modified  refine- 
ments based  on  strain  energy  and  strain  energy  density  shown 
in  Figure  5.2  (g)  and  (h)  .  As  can  be  seen  in  Table  I,  such  a 
modification  did  indeed  significantly  reduce  the  errors  in 
the  maximum  displacement  and  the  integral,  but  it  even 
further   improved   the    accuracy    for   the    maximum   slope. 

One  might  conclude  from  Table  I  that  the  best  overall 
model  was  obtained  using  the  graded  mesh,  and  that  since  it 
is  easier  to  obtain,  it  should  be  deemed  the  optimal  grid. 
But  this  particular  grading  was  caosen  for  precisely  that 
reason  and  was  presented  only  as  a  means  cf  comparison.  In 
practice,  the  selection  of  a  grading  ratio  is  somewhat  arbi- 
trary  and   making   an    adequate   choice    nay    be   difficult. 

There  is  justifiable  confusion  as  to  which  refinement 
produced  the  "best"  solution  accuracy  for  this  problem  and 
i-    raises   perhaps   the  most    important    issue    in   the    subject   of 
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grid  refinement.  Before  any  optimization  process  can  be 
pursued,  the  optimization  goals  must  be  explicitly  defined. 
Clearly,  as  is  the  case  in  this  problem,  the  designation  of 
the  optimum  grid  would  depend  heavily  upon  which  of  the 
three    solution    resultants   is   most    critical   to   the   analysis. 

B.       TAPEBED    BAR    PK03LEH 

The  linearly  tapered  bar  under  axial  loading  has 
received  considerable  attention  and  was  one  of  the  early 
problems  for  which  analytical  grid  optimization  was 
employed.  Consider  a  tapered  elastic  bar  of  length  L  and 
modulus  E,  fixed  at  one  end,  with  an  axial  load  P  applied  at 
the      other,  for      which        the      axial        displacement      u{x), 

(0  <  x  <  L)  ,  is  desired.  The  cross-sectional  area  varies 
linearly    from  AQ    at    the   fixed   end   to    At    at    the   tip,    as    shown 
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Figure  5.3        Linearly   Tapered   3ar    Under    Axial   Loading. 

in   Figure   5.3.         The   analytical      solution    and    finite   element 
approach   are  presented  in   Appendix   3. 
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One  of  the  significant  features  of  the  tapered  bar 
problem  is  that  the  maximum  stress  can  be  very  difficult:  to 
model  accurately,  and  it  is  for  precisely  these  problems 
exhibiting  large  strain  gradients  that  grid  optimization 
becomes      most      beneficial.  Interestingly,        the      stresses 

obtained  at  the  element  midpoints  are  exact  for  this 
problem,  and  the  difficulty  arises  from  the  inability  of  the 
constant  slope  shape  functions  to  model  the  maximum  stress 
occurring      at  the      boundary.  In      examining   this      problem, 

Prager  [Bef.  20]  demonstrated  analytically  that  when  each 
element  has  the  same  strain  energy  content,  the  relative 
error  in  displacement  is  identical  for  all  the  nodes. 
However,  this  phenomenon  appears  peculiar  to  this  problem 
and  the  author  does  not  subscribe  to  such  a  measure  of  an 
optimum   grid.  Judging   the      effectiveness   of     a    particular 

refinement  based  upon  the  deviation  or  the  mean  value  of  the 
pointwise  errors  generally  tends  to  be  unfavorable  to  opti- 
mization procedures  since  they  almost  always  introduce  many 
more  nodes  in  those  regions  where  the  response  is  most 
difficult  tc  model*  Hence,  an  improved  solution  may  have  a 
larger   mean    value  of   the   pointwise   errors   [Ref-    3  ]- 

The  criteria  employed  in  the  refinement  of  the  tapered 
bar  model  are  identical  to  those  used  in  the  cable  problem 
and  their  effects  are  shown  in  Figure  5.4  (a-e)  .  Two  excep- 
tions are  that  now  the  displacement  and  strain  energy 
criteria  produce  identical  refinements,  and  the  graded  model 
chosen  as  the  best  overall  is  now  based  on  a  grading  ratio 
of  1.4,  producing  a  more  drastic  refinement  than  that  of  1.2 
for  the  cable.  This  further  demonstrates  the  difficulty 
involved  in  obtaining  adequate  element  grids  on  an 
"a-priori"    basis. 

The  solution  results  are  presented  in  Table  II  and  the 
most  readily  apparent  observation  for  this  problem  is  the 
large   errors     in   the      maximum    slope,         which   would      severely 
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(a)  Uniform 


(b)  Graded  1.4 


(c)   u  and  U 


:d)  u 


(e)   SED 


(f)  SED   , 
mod 

Figure  5.4   Tapered  3 


ar  Sefinements. 


46 


TABLE    II 
Tapered    Bar   Problem   Solution  Results 

Problem  Parameters: 


L 

100 

in 

A 

10.5 

in2 

A 

0.5 

in2 

2 

=    10x1 06 

DSl 

P 

=    10x103 

lb 

Variation  Percentage    Relative   Error 

Refinement  fu 

Criterion 


Uniform 


j  Graded  (1.4) 

I  u    :   a 


I 

|  SED 


u (max) 

u'  (max) 

J0udx 

-3.80 

-37.5 

0.68 

-0.78 

-4.  1 

0.  14 

-0.85 

-  10.6 

0.  15 

-1.81 

-7.7 

0.33 

-6.54 

-3.6 

1.  18 

SED    (mod)  -1.99  -3.6  0.36 


underestimate  the  maximum  stress.  These  results  are  based  on 
quadratic  extrapolation  of  the  exaot  slopes  at  the  element 
midpoints,  since  the  linear  shape  functions  would  produce 
even  poorer  estimations  of  the  maximum  slope.  As  before,  the 
more  extreme  refinement  based  on  the  strain  energy  density 
variation  provides  the  most  accurate  estimation  of  the 
maximum  slope,  but  with  the  accompanying  degradation  in 
estimates  for  maximum  displacement  and  the  integral  of 
displacement.  Again,  the  large  errors  in  these  values  may 
be  significantly  reduced  by  employing  a  grading  scheme  to 
restrict  the  size  of  the  larger  elements  as  shewn  in  Figure 
5.4  (f)  .  Unlike  the  previous  problem,  such  a  modification 
has  no  effect  on  the  estimate  of  maximum  slope  because  of 
the  extrapolation  of  the  element  midpoint  slopes,  which  are 
exact   regardless   of    the   element   modal. 
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A  different  version  of  the  tapered  bar  problem,  for 
which  the  displacement  and  strain  energy  criteria  will  not 
produce      identical      refinements,  involves      replacing     the 

concentrated  tip  load  P  with  a  linearly  varying  axially 
distributed  load  q(x),  specified  by  the  values  at  the  fixed 
end  g^,  and  the  tip  qt  .  The  problem  may  be  further  modified 
by  reversing  the  bar  such  that  the  maximum  slope  occurs  at 
the    fixed   end,       while  the   maximum    displacement   occurs   at   the 
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Figure   5.5        Reversed  Tapered   Bar   with   Distributed  Load, 

free  end  as  shown  in  Figure  5.5.  The  case  of  the  linearly 
varying  distributed  load  is  included  in  the  formulation  in 
Appendix   B. 

This  problem  was  solved  for  a  uniformly  distributed  load 
using  the  same  procedure  as  in  the  previous  two  problems. 
The  refinement  models  are  presented  in  Figure  5.6  and  the 
solution  results  in  Table  III.  The  observations  are  consis- 
tent with  those  made  in  the  previous  problems,  but  now  one 
would    likely  agree      that    the  refinement    based      on    the    strain 
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TABLE    III 
Reversed    Tapered    Bar  Solution   Results 


Problem  Parameters:        L   =  100   in 

A    =  0.5   in 2 


Variation 

Refinement 

Criterion 

Uniform 

u 

u» 

0 

SED 

SED    (mod) 

A  =  10. 5  in* 
E  =  10x106  psi 
g  =    100  lb/in 

Percentage  Relative  Error 


u (max) 

u*  (max) 

J0udx 

-5.5 

-39.4 

-7.  1 

-2.0 

-7.5 

-2.6 

-2.7 

-3.  1 

-3.4 

-3.7 

-5.9 

-4.7 

-11.9 

-3.4 

-15.3 

-3.1      -3.4      -4.0 


energy  density  would  represent  an  optimal  grid,  provided 
that  modifications  are  introduced  to  prevent  any  elements 
from    growing  axcessively   large. 

C.       GUIDELINES    FOR    ONE-DIMENSIONAL    GRID    OPTIMIZATION 

The  most  important  lesson  to  be  learned  from  this  one- 
dimensional  study  is  that  the  grid  optimization  procedure  is 
necessarily  dictated  by  the  optimization  goal,  or  the  under- 
lying purpose  for  performing  the  finite  element  analysis.  No 
element  grid  oan  possibly  provide  optimum  accuracy  for  every 
solution  resultant  cf  interest.  In  solving  these  simple 
problems,  a  balance  has  been  sought  for  achieving  adequate 
accuracy        for      three        of      the        mora      important        solution 
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resultants,  with  emphasis  on  the  maximum  value  of  the  deriv- 
ative of  the  dependent  variable,  which  more  often  is  net 
only  the  mest  important  part  of  the  solution  but.  also  the 
most  difficult  to  obtain  accurately  in  finite  element 
analysis. 

The  important  grid  optimization  techniques  of  intro- 
ducing mere  degrees-cf-freedom  by  subdividing  the  elements 
or  increasing  their  polynomial  order  have  been  intentionally 
omitted  in  favor  of  the  optimization  strategy  of  seeking 
maximum  solution  accuracy  for  a  specified  number  of 
degr ees-of-f reedom  using  linear  elements.  This  is  because 
such  a  procedure  is  not  so  straight  forward  in  two- 
dimensional  problems  where  the  number  of  degrees-cf-f reedom 
are  dependent  on  some  geometric  considerations,  which  do  not 
appear  in  problems  of  one-dimension.  Based  on  this  choice  of 
optimization  strategy,  it  appears  the  strain  energy  density 
variation  provides  the  most  useful  criterion  for  the  adap- 
tive refinement  of  the  initial  grid.  Yet  all  three  problems 
demonstrated  some  pathological  results  that  can  arise  when 
the  elements  are  permitted  to  grow  excessively  large  in  the 
regions  where  the  strain  energy  density  varies  the  least.  In 
applying  a  scheme  tc  restrict  the  size  of  the  largest 
elements,  no  mention  has  been  made  of  how  tc  determine  when 
an  element  is  excessively  large.  It  has  become  the  experi- 
ence of  the  author  that  any  element  representing  over  half 
of  the  domain  should  probably  be  considered  too  large,  and 
measures   should    be    employed    to    restrict   its   size. 

It  would  appear,  at  least  for  these  classes  of  problems, 
that  this  difficulty  of  decreasing  accuracy  of  a  particular 
solution  parameter  for  successive  refinements  can  be  ignored 
by  merely  accepting  the  largest  value  among  the  cycles  as 
the  most  accurate  solution  result.  For  example,  it  was 
demonstrated  that  the  refinement  based  on  strain  energy 
density    provided   significant  improvement   in    the    accuracy    for 
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the  maximum  slope  but  underestimated  the  maximum  displace- 
ment even  mora  than  the  initial  uniform  grid.  Assuming  that 
the  linear  element  model  always  underestimates  such  maxima, 
the  maximum  slope  for  refined  grid  and  the  maximum  displace- 
ment for  the  unrefined  grid  could  be  accepted  as  the  op-imal 
results  of  the  analysis.  The  fallacies  of  such  a  procedure 
are  that,  first,  the  refinement  may  not  represent  the 
optimal  grid  as  it  has  been  defined  and,  second,  for  self- 
adaptive  finite  element  codes  the  user  is  provided  with  the 
"optimum  grid"  of  the  final  cycle  and  the  solution  results 
thereof. 

Melosh  and  Marcal  [Ref.  21]  have  proposed  an  alternative 
use  of  the  refinement  criterion  based  on  strain  energy 
density  variation  which  avoids  the  problem  of  excessive 
element  growth  altogether.  Beginning  with  a  reasonably 
coarse  unifora  grid,  those  elements  with  the  greatest  strain 
energy  density  variation  are  selectively  refined  by  either 
subdividing  them  or  increasing  their  polynomial  order  with 
the  introduction  of  additional  degrees-of-f reedom.  While 
such  a  procedure  does  not  equi-distribute  the  element  strain 
energy  variations,  it  can  reduce  them  all  to  some  prespeci- 
fied  tolerance,  such  as  a  percentage  of  the  average  element 
variations  for  the  initial  analysis.  Because  this  procedure 
is  particularly  attractive  for  grid  refinement  in  problems 
of  higher  dimensions,  it  will  be  employed  extensively  for 
the  study  of  grid  optimization  for  two-dimensional  problems 
in   the  next   chapter. 
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VI.     APPLICATION    TO    TWO-DIMENSIONAL    PROBLEMS 

Since  investigators  began  working  in  the  field  of  finite 
element  grid  optimization  in  the  early  1970's,  nearly  all  of 
the  effort  has  been  devoted  to  the  development  of  a  system- 
atic procedure  for  obtaining  optimal  grids  for  two- 
dimensional  problems  of  elasticity.  Even  today  there  are 
several  competing  approaches  to  this  problem  and  no  partic- 
ular one  has  yet  been  overwhelmingly  accepted  as  the 
preferred  method  of  grid  optimization.  While  it  is  the  two- 
dimensional  problem  for  which  most  of  these  techniques  have 
been  developed,  their  application  to  such  can  be  much  more 
difficult  than  for  the  one-dimensional  case.  Almost  invari- 
ably when  performing  grid  refinement  on  two-dimensional 
domains,  the  analyst  is  confronted  with  the  problems  of 
maintaining  int  ere  lenient  compatibility  and  preventing  severe 
element   distortion. 

In  selecting  an  appropriate  two-dimensional  problem  for 
the  application  of  seme  grid  optimization  techniques  and  a 
comparison  of  their  effectiveness,  it  is  desirable  that  the 
test    case   possess  the   following    properties: 

•  the  analytical  solution  should  exist  in  order  to 
perform    reliable    error    analysis 

•  the  solution  should  exhibit  sufficiently  large 
gradients  to  provide  a  meaningful  measure  of  the 
refinement    effectiveness 

•  the  idealization  should  have  one  degree-of-f r eedom  per 
node  and  possess  simple  boundary  conditions  to 
minimize  the  computational  effort  involved  in 
repeated  solutions 
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There  are  few  problems  that  meat  these  criteria,  but 
Saint-Venant  torsion  of  a  non-ciroular  section  provides  a 
good    test   case. 

A.       PROBLEM    DESCRIPTION 

Consider  a  solid  circular  shaft  of  radius  "a"  made  from 
isotropic  material  of  shear  modulus  G  and  having  a  circular 
groove,      or   keyway,       of   radius      b        along   a   generator   of   the 


Figure  6.1        Cross-section    of   Shaft   with  Keyway. 

shaft.    The   shaft  cross-section   is   shown    in    Figure    6.1.         The 
shaft   is   subjected      to  an   applied   torque        T      which    produces 
an   angle      of  twist    per   unit      length    9.        The    problem      may  be 
solved   by    finding      the  Prandtl   torsional    stress      function   $ 
which   satisfies    the    governing    differential   equation: 


.     2    +    .    2    +    2    =    0 

dX  3y 


(Eqn  .    6.1) 
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subject  to  the  Dirichlet  condition  that  ♦  =  0  on  the 
section  boundary.  The  torsional  stress  function  is  defined 
such  that  the  shear  stress  x  at  any  point  on  the  domain  may 
be    expressed  as: 

t     =G8     [(  §|,2+   (|i  ,2p  (Eqri.    6.2, 

For  this  formulation,  the  angle  of  twist  9  is  prescribed, 
rather  than  the  applied  torque  T.  The  torque  is  is  calcu- 
lated  from   the    area    integral: 

T    =    2G9  /     ♦    dA  (Eqn.    6.3) 


The  analytical  solutions  of  Equations  6.1  and  6.2  are 
derived      by        Sokoinikoff      [Sef.    22:    pp.    141-143]        and      are 

presented  in  Appendix  C  along  with  the  evaluation  of 
Equation  6.3  and  a  prescribed  finite  element  formulation. 
For  this  problem,  the  three  solution  resultants  of  interest 
for    the   grid  optimization   study   are: 


(1)  maximum  value  cf  the  dependent  variable,  or  torsion 
function      V^rnax; 

(2)  maximum  value  cf  the  gradient  of  the  dependent  vari- 
able (a  quantity  proportional  to  maximum  shearing 
stress       Tmax)  ; 

(3)  the  area  integral  of  the  dependent  variable  ever  the 
domain  (a  quantity  proportional  to  the  applied 
torque  T) . 

These  quantities  -  the  dependent  variable,  its  gradient,  and 
an  integral  thereof  -  are  selected  as  representative  of 
entities  whose  error  one  might  wish  to  minimize  in  a  finite 
element   analysis. 
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B.       COMPUTER   IMPLEMENTATION 

As  can  be  seen  in  Figure  6.1,  the  domain  of  this  problem 
is  symmetric  about  the  x-axis,  therefore  the  finite  element 
solution  need  only  be  obtained  for  the  upper  half  of  the 
domain.  For  all  of  the  solutions  presented  herein  the 
problem  geometry  is  defined  by  assigning  the  dimensionless 
ratio,  b/a  =3.4,  and  an  acceptable  upper  limit  on  the  anal- 
ysis cost  was  arbitrarily  chosen  to  be  that  corresponding  to 
approximately  one  hundred  nodal  points.  The  computation  and 
assembly  of  the  finite  element  matrices  and  solution  of  the 
resulting  system  of  equations  was  performed  using  the  steady 
state  heat  conduction  operations  of  CAL-NPS  [  Ref -  23].  This 
group  of  subroutines  comprises  an  efficient  finite  element 
code  for  solving  Poisson's  eguation  in  two  or  three  dimen- 
sions and  has  the  additional  advantage  of  permitting 
variable-noded    isoparametric  elements. 

Since  there  was  no  readily  available  interactive  prepro- 
cessor which  lent  itself  well  to  adaptive  mesh  refinement, 
the  author  had  no  choice  but  to  create  his  own.  Since  the 
problem  domain  is  simply  connected,  the  automatic  mesh 
generation  was  performed  employing  inverse  mapping  of  a 
single  cubic  isoparametric  element  of  the  serendipity  family 
onto      the      problem    dcmain      [Ref.    24:    pp.    228-229].  Mapped 

boundary  nodes  were  repositioned  to  conform  to  the  actual 
domain  boundary  and  additional  nodes  generated  during  the 
refinement    process    were   mapped    using    the   same    procedure. 

Since  the  finite  element  code  selected  for  this  investi- 
gation provided  output  only  for  the  nodal  values  of  the 
dependent  variable,  it  was  coupled  to  the  author's  postpro- 
cessor. Such  a  postprocessor  is  necessary  in  the  optimiza- 
tion process  for  computing  nodal  values  of  shear  stresses 
and  strain  energy  density,  element  contributions  to  torque 
and    total   strain   energy,    and   exact   results    from    theory. 
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C.       ASYMPTOTIC    ERROR    ANALYSIS    FOR    UNIFORM    REFINEMENT 

The  concept  of  asymptotic  convergence  rate  for  uniformly 
refined  grids  was  presented  in  Chapter  3.  When  the  number  of 
uniformly  distributed  degrees-of -freedom  is  sufficiently 
large,  the  log-log  plot  of  the  relative  energy  error  versus 
the  number  of  degrees-of -freedom  is  approximately  linear  in 
the  asymptotic  range.  The  slope  of  this  line  represents  the 
asymptotic   rate    of    convergence   in   energy. 

It  so  happens  that  relative  error  in  the  torque  T  of 
this  problem  is  equal  to  the  relative  energy  error  and 
therefore  exhibits  this  linear  asymptotic  behavior  on  the 
log-log  plot  against  the  number  of  uniformly  distributed 
degrees-of- freedom.  Fortunately,  the  other  two  solution 
resultants  of  interest  behave  similarly.  This  will  prove 
very  beneficial  in  performing  the  error  analysis  for  this 
two-dimensional  study  for  two  reasons.  First,  because  it  is 
unnecessarily  difficult  to  construct  an  optimal  grid  with 
the  same  number  of  degrees- of- freedom  as  a  uniform  grid,  the 
linear  behavior  of  the  solution  resultants  in  the  asymptotic 
range  on  the  log-log  scale  permits  interpolation  for  any 
number  of  degrees-of -freedom.  Then  the  solution  results  for 
a  uniform  grid  of  the  identical  number  of  degrees-of-f reedom 
provides  a  reference  for  comparison  to  determine  the  effec- 
tiveness of  the  optimization  technique.  Secondly,  if  the 
convergence  rate  of  a  particular  solution  resultant  is 
extremely  slow,  as  is  often  the  case  for  maximum  stress,  it 
becomes  difficult  to  gain  an  appreciation  of  the  true  effec- 
tiveness of  the  optimization.  For  example,  an  order  of 
magnitude  reduction  in  the  relative  solution  error  may 
require  an  order  of  magnitude  increase  in  the  number  of 
degrees-of- freedom  using  uniform  refinement,  but  relatively 
few  additional  degrees-of -freedom  using  an  optimization 
technique.    Therefore,    it    will   be    enlightening   to    extrapolate 
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(a)    Linear   Element    Grid 


(b)     Quadratic    Element    Grid 


Figure    6.2        Uniform  Linear   and    Quadratic   Element    Grids 
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the  relative  error  versus  degrees-of- freedom  curve  to  obtain 
a  rough  approximation  of  the  number  of  degrees-of-f reedom 
necessary  to  obtain  solution  accuracies  similar  to  the 
optimal  grid,  but  using  successively  finer  uniform  grids.  Of 
course,  this  is  only  an  estimation  and  ignores  such  reali- 
ties as  numerical  ill-conditioning  and  computer  rcund-off 
error. 

A  uniform  grid  is  one  for  which  all  of  the  elements  are 
of  the  same  size  h  and  the  same  polynomial  order  p.  Clearly, 
it  is  impossible  to  obtain  such  a  grid  for  this  particular 
domain  using  isoparametric  mapping,  but  a  nearly  uniform 
grid  may  be  constructed  in  which  the  elements  are  of  approx- 
imately the  same  size.  Such  unifora  grids  are  shown  for  the 
cases  of  linear  quadrilateral  elements  and  quadratic  seren- 
dipity elements  (Fig.  6.2).  For  this  geometry,  the  uniform 
grid  is  not  uniquely  defined  for  a  specified  number  of 
elements.  This  is  because,  in  performing  isoparametric 
mapping,  there  must  be  specified  four  keynodes  on  the  actual 
domain  boundary  to  correspond  to  the  four  corner  nodes  of 
the      parent      square.        Since      this      domain      has      only      three 
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Figure   6.3        Keynode  Placement   for    Isoparametric   flapping. 


59 


vertices,  the  placement  of  the  fourth  keynode  is  at  the 
discretion  cf  ths  analyst  (Fig.  6.3),  and  can  have  a  notice- 
able   effect  on    the    solution    results. 

The  asymptotic  error  analysis  was  performed  for  the 
three  solution  resultants  of  interest  using  uniform  grids  of 
linear  and  quadratic  elements.  The  results  are  presented  in 
Figure  6.4.  All  of  the  solution  resultants  behaved  as 
predicted  with  the  exception  of  the  maximum  torsion  function 
value   using   linear    elements,      y/m  (1)  .  It   appears   that   the 

accuracy  of  this  particular  parameter  is  very  strongly 
dependent  upon  the  keynode  placement.  The  curve  constructed 
in  Figure  6.4  represents  an  average  for  several  keynode 
positions. 

While  Figure  6 . 4  is  intended  primarily  to  serve  as  a 
reference  tool  for  future  analyses,  it  provides  some  inter- 
esting  irformation: 

("!)  For  the  cases  of  maximum  torsion  function  value  and 
applied  torque  (and  energy),  the  asymptotic  rate  of 
convergence  using  quadratic  elements  is  more  than 
twice   that    for    linear    elements. 

(2)  While  the  error  in  torque  for  the  quadratic  case  is 
always  smaller  than  that  for  the  linear  case,  the 
linear  grid  may  provide  better  accuracy  for  the 
maximum  torsion  function  value  S^max  in  the  pre- 
asymptotic   range. 

(3)  3oth  accuracy  and  convergence  rate  in  the  maximum 
shear  stress  are  only  minutely  greater  for  the  quad- 
ratic  element    grid   than    for   the    linear   one. 

However,  for  this  last  observation,  the  point  must  be  made 
that  for  the  linear  element  grid,  the  maximum  shear  stress 
was  obtained  by  quadratic  interpolation  rather  than  from  the 
linear   shape  functions.      While   this    will   greatly    improve   the 
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accuracy  cf  the  maximum  shear  stress  approximation,  it  will 
have  no  effect  on  its  rate  of  convergence.  Therefore,  if 
obtaining  an  optimal  estimate  of  the  maximum  shear  was  the 
purpose  of  the  analysis,  there  is  much  to  be  said  on  behalf 
of  linear  elements  besides  their  computational  efficiency. 
Of  course,  this  observation  is  based  on  uniform  grid  refine- 
ment, which  would  rarely  compete  favorably  with  the  optimi- 
zation  techniques  to    be   examined. 

The   reason    that    the   rate   of  convergence   in   maximum    shear 
stress   is   so  poor   using      uniform   refinement   for   this   problem 


TCA)=0 


TC8)^TT 


max 


Figure  6.5        Stress  Distribution   on   Shaft  Keyway. 

can  be  seen  in  Figure  6.5.  The  shear  stress  varies  greatly 
over  a  short  distance,  by  increasing  from  zero  at  point  A  to 
its  maximum  value  at  point  B.  As  a  result,  there  exists  a 
region  of  excessively  large  strain  gradients  along  the 
keyway  which  severely  hinders  the  rate  of  convergence  when  a 
uniform  grid  is  employed.  If  the  keyway  radius  were  allowed 
to   approach   zero  producing    a   singularity   in   the    solution,    it 
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would  likely  be  necessary  to  employ  even  higher  crdsr 
elements  via  the  p-version  in  order  to  achieve  convergence 
using   uniform  refinement  [Ref.    1  ]. 

D.       PROBLEM    SOLUTION    WITH   GRID    OPTIMIZATION 

The  finite  element  solution  of  the  torsion  problem  will 
be  obtained  employing  the  following  grid  optimization  tech- 
niques  as   presented    in  Chapter   4: 

•  Contouring 

-  contours    of  the  torsion    function;    linear  elements 

-  contours    of   shear   stress;    linear   elements 

-  contours    of   strain    energy   density;    linear   elements 

•  Selective  Refinement 

-  h-version;    linear   elements 

-  h- version;    quadratic    elements 

-  p-version 

•  Subdcmain  Isolation 

-  linear  elements 

-  quadratic    elements 

•  Mesh   Grading 

-  linear  elements 

-  quadratic    elements 

1 .      Contouring 

The  original  finite  element  analysis  was  performed 
en  a  uniform  grid  of  98  linear  eieients,  78  nodes,  and  72 
degrees-of-freedom.  The  fiaite  element  solution  provided  the 
nodal  values  of  the  torsion  function  if  ,  from  which  the 
conventional  nodal  resultants  of  shear  stress  x  ,  and 
applied   torque      I      were   computed.    Based    upon    the    maximum   and 
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minimum  values  obtained  for  each  parameter,  along  with 
consideration  for  their  values  along  the  boundary,  the 
contours  to  be  used  for  nodal  placement  in  each  case  were 
selected.  The  number  of  contours  for  each  case  was  chosen  to 
maintain  approximately  the  same  number  of  degrees-cf-freedom 
as    for  the    initial    analysis. 

The  points  for  each  contour  value  selected  were 
obtained  by  linearly  interpolating  between  the  nodal  values 
of  each  parameter  obtained  from  the  initial  analysis.  The 
contours  were  constructed  by  smoothly  connecting  the  points 
by  hand.  The  element  layout  along  the  contours  posed  the 
most  formidable  problem  because  the  coarse-to-fine  tran- 
sition often  resulted  in  severe  elsaient  distortion,  and  it 
sometimes  became  necessary  to  degenerate  guadriiateral 
elements  into  triangles  when  the  transition  was  acute.  It 
was  decided  that  the  optimal  element  shapes  should  be 
preserved   along    the    contours   in   regions    of   highest    stress. 

The  contours  obtained  and  the  corresponding  grid  are 
presented   for  each    of  the   following    solution   resultants: 

•  torsion    function  (Fig.    6.6) 

•  shear   stress  (Fig.    6.7) 

•  strain   energy  density,    SED      (Fig.    6.8) 

The  resulting  grid  for  each  of  the  response  function 
contours  produces  smaller  elements  in  the  region  of  greatest 
stress  near  the  bottom  of  the  Jceyway  and  around  the 
periphery  of  the  shaft  where  the  stress  is  moderately  high. 
Consequently,  the  elements  near  the  center  where  the  stress 
is  zero  are  larger.  These,  of  course,  are  the  desired 
effects  for  an  optimization  criterion.  A  somewhat  unusual 
behavior  is  observed  at  the  point  of  intersection  of  the 
keyway  and  the  shaft  boundary  where  the  stress  is  also  zero. 
Apparently,  the  shear  stress  gradient  is  larger  than  the 
gradients   in  torsion    function    and  the   strain   energy    density, 
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(a)  Contours 


(b)  Corresponding  Grid 


Figure  6.6    Contouring  for  tas  Torsion  Function, 


6  5 


(a)  Contours 


(b)  Corresponding  Grid 


Figure  6.7   Contouring  for  the  Shear  Stress  Functi 


on. 
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(a)  Contours 


(b)  Corresponding  Grid 


Figure  6.8    Contouring  for  the  SED  Function, 
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resulting  in  smaller  elements  being  produced  in  that  region 
by  the  shear  stress  criterion  (Fig.  6.7)  than  by  the  other 
criteria.  While  all  of  the  grids  possess  to  some  extent  the 
desirable  features  cf  an  optimal  grid,  the  strain  energy 
density  function  produces  a  far  more  drastic  refinement 
towards  the  point  of  maximum  stress,  while  the  others  repre- 
sent more  moderate  refinements.  In  fact,  the  SED  contours 
are  sc  dense  around  the  keyway  that  the  coarse-to-fine 
element  transition  scheme  must  include  degenerate  quadrilat- 
erals to  avoid  violating  the  contours.  Note  also  that  the 
coarse-tc-f ine  transition  for  the  torsion  function  response 
is  fairly  smooth  whereas  the  transition  for  the  strain 
energy  density  and  shear  stress  refinement  tends  to  produce 
distorted  and  elongated  elements.  This  is  aggravated  oy  the 
fact  that,  unlike  the  torsion  function,  the  shear  stress  and 
strain   energy  density   are  not    monotonic   over   the    domain. 

The  solution  results  obtained  for  each  grii  are 
presented  in  the  upper  half  of  Table  IV.  At  first  glance, 
the  results  of  the  refinements  are  disappointing  in  compar- 
ison to  the  parenthetic  values  obtained  using  a  uniform 
grid,  while  all  three  criteria  produce  improved  accuracy  for 
the  maximum  shear  stress,  the  errors  for  the  maximum  torsion 
function  value  grow  extremely  large.  Recalling  the  observa- 
tions made  for  the  cne-dimensional  study,  it  would  follow 
that  such  behavior  is  probably  due  to  the  unusually  large 
elements  near  the  center  of  the  domain.  The  entries  in  the 
lower  half  of  Table  IV  reflect  the  drastic  improvement 
obtained  by  simply  introducing  a  few  additional  degrees-of- 
freedom  along  those  element  edges  which  grew  exceptionally 
long  during  the  optimization  process,  thus  increasing  their 
polynomial  order  from  one  to  two.  Not  only  did  this  modifi- 
cation reduce  the  large  errors  for  the  maximum  torsion 
function  estimation,  but  modest  improvements  were  also 
obtained   in    the    estimations    of    the   other   resultants    as    well. 
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Onca  again,  the  selection  of  the  optimum  grid  would  depend 
predominantly  on  the  optimization  goal  of  the  analysis,  but 
one  would  likely  agree  that  the  strain  energy  density  varia- 
tion along  with  some  modification  to  restrict  excessive 
element    growth    provides   the    superior    refinement    criterion. 

An  additional  word  of  caution  is  in  order  for  the 
contouring  techniques  for  grid  optimization.  Because  the 
problem  must  be  completely  redefined  from  scratch  after  the 
initial  analysis,  the  preprocessing  cost  can  become  enor- 
mous, especially  if  several  cycles  are  employed  to  obtain 
more  precise  contours  as  some  authors  suggest.  Unless  there 
is  available  an  interactive  automatic  mesh  generator  based 
on  this  technique,  such  as  the  one  described  in  Reference 
[Ref.  11  ]r  contouring  should  be  abandoned  in  favor  of  some 
more  easily  implemented  grid  optimization  techniques 
employing   similar   refinement  criteria. 

2-      Selective  Refinement 

The  simplest  way  to  avoid  the  problems  encountered 
in  the  contouring  techniques  is  to  perform  the  initial  anal- 
ysis on  a  reasonably  coarse  grid  and  then  to  selectively 
refine  these  elements  over  which  the  strain  energy  density 
varies  the  most.  The  critical  concern  then  arises  as  to  how 
coarse  the  initial  grid  should  be.  If  the  preprocessor 
employs  the  necessary  constraint  conditions  to  permit  the 
"fat  her-tc- four- sons"  element  subdivision  scheme  directly, 
or  if  hierarchical  refinement  is  employed,  then  the  initial 
grid  should  be  just  coarse  enough  to  adequately  define  the 
problem  and  to  limit  the  number  of  refinement  cycles  neces- 
sary. The  latter  becomes  even  less  of  a  concern  if  iterative 
solvers  are  employed.  If,  on  the  other  hand,  coarse-tc-f ine 
transition  schemes  are  used  to  implement  the  h-versicn  or 
only  lew  polynomial  order  elements  are  available  in  the 
p-version,      then   the    initial   grid      must    be    sufficiently   fine 
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so  as  net  to  restrict  severely  the  amount  of  refinement 
which  can  be  performed  in  any  givan  cycle.  Unfortunately, 
tha  conditions  under  which  this  investigation  was  conducted 
were    those   of  the   latter. 

a.      The    h- Version 

Selective  refinement  by  the  h-version  was 
performed  on  both  linear  and  quadratic  element  grids.  For 
the  linear  case,  the  initial  analysis  was  performed  on  a 
uniform  grid  of  55  elements,  72  nodes,  and  50  degrees-of- 
freedem.  The  initial  quadratic  analysis  employed  an  eight 
element  uniform  grid  of  37  nodes  and  20  degrees-of-freedom. 
The  reason  for  such  a  gr  eat.  disparity  in  the  number  of 
elements  for  the  initial  analyses  is  that  subdividing  a 
quadratic  element  introduces  many  more  degrees-of-freedom 
than  the  subdivision  of  a  linear  element.  Thesa  numbers  ware 
chosen  to  provide  approximately  tha  same  number  of  dagrees- 
of-freedem  for  the  optimum  grid  of  tha  final  cycle  for  each 
case.  The  initial  analysis  is  performed  and  those  elements 
over  which  the  strain  energy  density  variation  is  signifi- 
cantly greater  become  candidates  for  refinement.  The  refine- 
ment is  performed  by  subdividing  each  candidate  element  into 
four  new  ones  by  constructing  a  ooarse-to-fine  transition 
zone  of  "buffer"  elements  around  the  refined  region. 
Successive  analyses  and  selective  refinements  are  repeated 
until  the  maximum  element  strain  energy  density  variation  is 
approximately  that  of  the  remainder  of  tha  grid.  The  process 
is  further  improved  when  the  nodal  values  of  the  strain 
energy  density  are  used  to  indicate  the  general  direction  in 
which  the  refinement  is  to  proceed.  This  permits  multiple 
refinements  in  the  same  cycle,  thereby  reducing  the  number 
of  cycles  required  tc  arrive  at  tha  optimal  grid.  For  this 
problem,  the  linear  grid  required  two  refinement  cycles 
while      the      quadratic     grid      required      three      cycles.  The 
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(a)    Initial   3rid 


(b)     Cycle    1 


(c)     Cycle    2 


Figure   6.9      Sslective  Refinement    PuDceiure    -    Linear    Element: 
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selective  refinement  process  is  depicted  in  Figures  6.9  and 
6.10  for  the  linear  and  quadratic  element  grid, 
respectively. 

The  solution  results  for  each  selective  refine- 
ment cycle  are  presented  in  Tabls  V.  The  most  impressive 
observation  to  be  made  is  the  significant  improvement  in  the 
maximum  shear  stress  estimate  for  successive  cycles.  While 
there  is  also  modest  improvement  in  the  accuracy  of  the 
torque  estimate  for  successive  cycles,  when  compared  to  the 
estimate  obtained  from  the  uniform  grid  of  the  same  number 
of  degrees-of-f reedom  and  polynomial  order,  the  refinement 
estimate  cf  torque  is  slightly  poorsr.  This  is  because  addi- 
tional degrees -of- freedom  are  being  introduced  in  only  a 
small  region  of  the  domain  but  the  torque,  and  energy,  are 
global  quantities.  The  author  has  no  satisfactory  explana- 
tion why  the  estimate  for  the  maximum  torsion  function 
improves  for  successive  refinements  of  the  linear  grid  but 
not  for  the  quadratic  case.  Howevsr,  as  has  already  been 
mentioned,  this  particular  solutioQ  parameter  appears  very 
sensitive  tc  such  problem  variables  as  nodal  placements  and 
element  shapes;  hence,  its  behavior  is  difficult  to  predict, 
even  when  the  refinement  is  applied  to  regions  remote  from 
the  point  where  the  maximum  torsion  function  value  occurs, 
as      was    the     case      in     these   examples.  For      computational 

reasons,  it  is  desirable  to  restrict  the  number  of  refine- 
ment cycles  to  a  necessary  minimum.  In  this  example,  the 
quadratic  grid  required  an  additional  cycle  over  the  linear 
grid  but  this  is  because  it  is  necessary  to  perform  the 
initial  quadratic  analysis  using  far  fewer  degrees-of- 
freedom.  Therefore,  the  early  cycles  of  the  quadratic  anal- 
ysis   actually   represent   comparatively   smaller    problems. 
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b.      The    p-7ersion 

Before  continuing  to  ths  next  optimization  tech- 
nique, it  is  worthwhile  to  take  a  quick  look  at  selective 
refinement      employing  the      p-version.  Because    the      finite 

element  code  used  in  this  investigation  only  provided  for 
element  orders  of  one  and  two,  the  advantages  of  the  method 
cannot  be  fully  realized,  but  the  affects  of  a  single  cycle 
can    be   examined. 

Beginning  with  three  uniform  grids  with 
differing  numbers  of  linear  elements,  the  initial  analyses 
were  performed.  In  each  case,  the  elements  over  which  the 
strain  energy  density  varied  the  most  were  transformed  from 
4-noded  Lagrangian  elements  into  8-noded  serendipity 
elements  by  the  addition  of  midsids  nodes.  The  element 
grids  are  shown  in  Figure  6.11  and  the  asterisks  denote 
those  elements  for  which  the  polynomial  order  was  increased. 
In  this  example,  the  number  of  elements  to  be  refined  in 
each  case  was  chosen  so  as  to  achieve  approximately  the  same 
number  of   degrees-of-freedom  after  a    single   cycle. 

The  solution  results  are  shown  in  Table  71. 
Significant  improvements  in  the  estimate  of  ths  maximum 
shear  stress  were  achieved  for  each  case.  An  improvement  in 
the  estimated  torque  was  also  realized  for  all  three  cases, 
but  was  more  noticeable  when  the  number  cf  refined  elements 
was  larger.  This  is  because  quadratic  elements  are  far 
superior  to  linear  elements  in  ;as  modeling  of  integral 
quantities,  as  observed  in  Figure  6.4.  Somewhat  disturbing 
is  the  increased  error  in  the  estimate  of  the  maximum 
torsion  function  value  observed  in  two  of  the  three  refine- 
ments even  though  the  elements  in  the  vicinity  of  its  occur-* 
rence  were  net  affected.  Again,  this  is  likely  attributable 
to  the  unusual  behavior  patterns  of  this  quantity  already 
discussed. 
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(a)    76   Degrees-of-Fre edcm   -   23    Rafined   Elements 


(b)    76   Degrees -of-Fre eaom    -    9    Refined    ^ieaer.is 


(c)    75   Degrees-of-Freedoni    -    4    Refined    Elements 


Figure   6.11        Selective   Refinement    Employing   p-Version. 
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In  order  to  perform  additional  cyclas  of  the 
p-version,  it  would  be  necessary  to  alter  the  refinement 
criterion  slightly.  Because  the  element  sizes  do  not  change 
for  successive  cycles,  the  need  for  rafineraent  would  neces- 
sarily be  based  on  strain  energy  density  variation  between 
nodes   rather  than  over  the    elements. 

Selective  refinement  employing  the  p-version  is 
most  efficiently  implemented  hierarchically,  in  which  case 
it  acquires  some  attractive  computational  advantages.  It  is 
unfortunate  that  time  did  not  permit  further  investigation 
here,    but   the  need    for  future    research   is    evident. 

3 .      Subdomain  Isolation 

The  refinement  criterion  and  initial  procedures  in 
employing  subdomain  isolation  are  identical  to  those  used  in 
selective  refinement.  After  the  candidate  elements  for 
refinement  are  identified,  they  ara  completely  isolated  from 
the  remainder  of  the  domain  and  solved  as  smaller  subdomain 
problems.  The  advantages  of  the  technique  are  twofold.  By 
isolating  the  elements  to  be  refined  the  solution  is  not 
repeated  in  each  cycle  for  those  alaments  for  which  the 
initial  analysis  is  assumed  adequate.  Furthermore,  by  elimi- 
nating most  of  the  degrees-of -freedom  over  the  entire 
domain,  the  subsequent  refinement  of  the  isolated  region  can 
be   much   greater    than    would    otherwise    be    practical. 

As  before,  the  technique  was  applied  to  both  linear 
and  quadratic  uniform  element  grids.  Those  elements  of  the 
initial  analysis  over  which  the  strain  energy  density  varia- 
tion was  exceptionally  large  were  isolated  to  comprise  the 
subdomain  in  each  case.  There  were  three  such  elements  of 
the  initial  linear  grid  and  two  for  the  quadratic  grid.  Each 
subdomain  was  uniformly  refined  to  achieve  approximately  the 
same  number  of  degrees-cf- freedom  as  the  initial  analyses. 
The  process  is  depicted  in  Figure  5.12  for  the  linear  grid 
and    Figure   6.13    for    the   quadratic   case. 
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(a)    Initial   3rid 


(b)     Refineoisnt    of    Isolatsd    Subdotnain 


Figure    6.12        Subdomain    Isolation    -    Linear   Elements, 


8  0 


(a)    Initial   Srid 


(b)     Refinement    of    Isolated  Subdomain 


Figure   6.13        Subdomain   Isolation    -    Quadratic    Elements. 
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In  performing  sub  domain  isolation  the  governing 
equations  remain  the  same,  while  only  the  domain  and  the 
boundary      conditions     are      altered.  When     the      subdomain 

boundary  has  nodes  common  to  the  initial  grid,  then  -he 
boundary  values  for  those  nodes  are  simply  the  solution 
values  obtained  from  the  initial  analysis.  The  boundary 
values  arising  from  the  introduction  of  new  boundary  nodes 
during  the  refinement  process  must,  be  generated  by  interpo- 
lation of  the  solution  results  of  the  initial  analysis.  One 
of  the  options  for  an  interpolation  scheme  is  simply  to  use 
the  shape  functions  of  the  unrefined  elements.  This  may  not 
be  desirable  in  the  case  of  linear  elements,  so  a  higher 
order  interpolation  may  be  employed.  In  this  example  a  third 
order  Lagrangian  polynomial  was  convenient  for  the  linear 
case  since  there  are  four  nodes  from  the  initial  analysis 
along  the  right-hand  boundary  of  the  subdomain.  Since  there 
are  only  two  such  nodes  on  the  upper  boundary,  it  is  neces- 
sary to  "borrow"  some  adjacent  nodal  values  from  the 
discarded  portion  of  the  domain  in  the  generation  of  new 
boundary    values    using  higher   order  interpolation. 

The  solution  results  presented  in  Table  VII  are 
quite  remarkable.  In  a  single  cycle,  the  solution  accuracy 
for  the  maximum  shear  stress  has  increased  by  a  full  order 
of  magnitude.  No  other  optimization  technique  examined  in 
this    investigation      produced  such   improvements.  Note   that 

the  higher  order  polynomial  interpolation  for  the  boundary 
values  did  improve  the  solution  results  for  the  linear  case. 
One  of  the  disadvantages  of  this  technique  is  that  the 
refinement  can  produce  no  improvement  in  the  estimation  of 
local   quantities   outside   the   subdomain.       As   in   this    example, 


the 


estimation 


the 


maximum      torsion      function      value 


obtained  from  the  initial  analysis  must  be  accepted  as  the 
optimum  since  it  occurs  outside  the  subdomain.  Furthermore, 
since   its    value    is    predominantly      affected    by   refinements   in 
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regions  of  high  gradients,  it  is  doubtful  that  isolating  a 
new  subdomain  using  the  elements  adjacent  to  its  point  of 
occurrence  would  noticeably  improve  its  accuracy.  However, 
since  the  torque  is  a  globally  computed  quantity,  refinement 
will  improve  the  accuracy  of  its  contribution  from  the 
subdomain  resulting  in  improvement  of  its  accuracy  overall 
as  observed  in  Table  VII.  It  is  this  strictly  local  nature 
of  the  subdomain  isolation  technique  which  restricts  its 
applicability.  But  if  the  optimization  goals  are  well 
defined  and  it  is  understood  under  which  conditions  and  for 
what  parameters  it  is  effective,  it  can  be  an  extremely 
powerful    grid  optimization    technique. 

4  •      &§s  h  3rading 

While  mesh  grading  is  nearly  always  performed  on  an 
"a-priori"  basis,  it  may  also  be  employed  adapt ively  to 
provide  a  simple  grid  optimization  technique.  After  an 
initial  solution  has  been  obtained,  the  analysis  may  be 
repeated  using  various  combinations  of  grading  ratios  in 
order  to  achieve  a  more  uniform  distribution  of  the  element 
strain  energy  density  variations.  Here  the  grading  ratio 
refers  to  the  constant  ratio  of  adjacent  element  lengths 
along  a  boundary  of  the  domain  to  which  grading  is  applied- 
There  are  several  drawbacks  to  the  technique,  the  first 
being  that  a  good  combination  of  grading  ratios  may  be 
difficult  to  obtain  in  a  reasonable  number  of  cycles.  The 
other  difficulty  is  that  if  smaller  elements  are  desired  in 
more  than  one  region  the  domain  must  be  refined  and 
constructed    by   subregions. 

Unfortunately,  the  domain  of  this  problem  is  not 
well  suited  for  mesh  grading  since  it  possesses  no  favorable 
geometric  symmetry,  Hence,  the  resulting  element  elongation 
and  distortion  would  become  severe  for  larger  grading 
ratios.      For  simplicity,    the   nodal   placements    will    be    biased 
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only  along  tha  two  domain  boundaries  adjacent  to  the  point 
of  maximum  stress  and  the  same  grading  ratio  will  be  applied 
for  each.  This  will  result  in  small  slements  near  the  bottom 
of  the  keyway  and  large  elements  along  the  periphery  of  the 
shaft.  While  this  is  not  the  most  desirable  grid  topology, 
it  will  produce  a  more  uniform  distribution  of  the  element 
SED    variations. 

The  technique  was  applied  to  both  linear  and  quad- 
ratic element  grids  starting  with  a  uniform  mesh  and  succes- 
sively increasing  the  grading  ratio  until  the  elements  along 
the  shaft  periphery  exhibited  S£D  variations  as  large  as 
those  for  the  elements  along  the  keyway.  In  both  cases  this 
condition  occurred  beyond  the  point  where  excessive  element 
elongation  would  be  expected  to  produce  numerical  inaccura- 
cies. Graded  meshes  for  selected  grading  ratios  r,  are  shown 
in  Figure  6.14  for  linear  elements  and  Figure  6.15  for  quad- 
ratic elements.  The  solution  rasults  are  presented  in 
Taole  VIII.  As  can  be  observed,  the  maximum  shear  stress 
estimate  improves  for  each  successive  increase  in  the 
grading  ratio.  However,  the  cost  of  such  improvement  is  the 
accompanying  degradation  in  the  sstimate  of  the  maximum 
torsion  function  value.  This  is  to  be  expected  since  the  two 
maxima  occur  at  different  locations  in  the  domain  and  there- 
fore decreasing  the  size  of  the  elsments  in  the  vicinity  of 
one  will  necessarily  increase  the  element  size  near  the 
other.  Note  that  this  degradation  is  not  nearly  so  severe 
for  the  quadratic  elements,  and  the  accuracy  actually 
improves  for  a  low  valus  of  the  grading  ratio.  This  is 
because  the  higher  order  interpolation  can  better  accomodate 
larger  elements.  For  both  linear  and  quadratic  element  grids 
the  optimal  accuracy  in  torque  estimation  occurs  for  the 
moderately   graded   meshes. 
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TABLE   VIII 

Hesh   Grading 

Solution    Results 

Perc 

entage   Relative   Error 

Grading                 Wy 
Ratio                       rmax 

Tmax    /G9             I    /    G9 

Linear  Element   Grids 

(    78    elsments,    93    nodes, 
72    degress-of- freedom   ) 

1.0                          0.060 

-6.06                      -1.77 

1.1                           0.161 

-2.63                      -1.56 

1.2                          0.339 

-1.03                      -1.77 

1.3                          0.679 

-0.484                    -2.20 

Quadratic   Element  Grid 

s    (    28    elements,    107   nodes, 
76    dagraes-of-freedom    ) 

1.0                        -0.0093 

-5.26                     -0.0116 

1.1                          0.0063 

-1.85                     -0.0064 

1.2                          0.0162 

-0.606                   -0.0091 

1.3                       -0.0246 

-0.413                   -0.0179 

I 


It  is  likely  that  some  of  the  error  is  attributable 
to  the  numerical  inaccuracies  dua  to  element  distortion. 
When  applying  a  grading  technique  the  analyst  should  seek  an 
equitable  balance  between  the  refinement  criterion  and  the 
grid    topology. 
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Finally,  since  the  grading  ratio  is  usually  applied 
to  the  nodal  separation  rather  than  the  element  edge 
lengths,  it  is  advisable  to  reposition  the  midside  nodes  of 
quadratic  elements  so  that  they  lie  near  the  center  of  the 
element  edges.  This  will  generally  improve  the  accuracy  of 
all  the  solution  parameters,  especially  if  the  grading  is 
somewhat    extreme. 

E.       GUIDELINES    FOB    T SO-DIMENSIONAL    5RID    OPTIMIZATION 

In  order  to  provide  some  guidelines  for  obtaining 
optimal  finite  element  solutions  for  two-dimensional  prob- 
lems it  is  helpful  tc  compare  the  solution  results  obtained 
for  this  problem  employing  the  optimization  techniques 
available.  Such  a  comparison  is  presented  in  Table  IX.  The 
upoer  portion  is  for  those  techniques  for  which  the  initial 
analysis  was  performed  using  linear  elements  and  the  lower 
portion  using  quadratic  elements.  Note  that  ail  of  the  grids 
employ  approximately  the  same  number  of  degrees-of- freedom, 
which  was  the  chosen  measure  of  analysis  cost  ic  this 
investigation. 

In  making  this  comparison  it  is  important  to  understand 
just  how  significant  a  change  in  error  actually  is.  If  the 
convergence  rate  of  a  solution  paramstar  is  very  slow,  even 
a  small  reduction  in  the  error  may  require  many  more 
degrees-of-  freedom.  For  this  reason,  the  numbers  in  paren- 
theses have  been  included  by  aach  arror  to  provide  a  rough 
approximation  of  the  number  of  dagrees-of-f reedom  required 
to  obtain  similar  accuracy  using  a  uniform  grid  of  the  sane 
number  of  degrees-of -freedom  and  elaments  of  the  same  poly- 
nomial order  as  tha  initial  analysis.  In  some  cases,  the 
analyses  were  not  actually  performed  but  instead  the  numbers 
in  parentheses  were  obtained  by  extrapolation  of  the  error 
versus  degrees-of-f reedom  curve  (Fig.  6.4),  and  ignoring 
round-off  and   ill-conditioning   effects. 
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The  first  observation  to  be  mais  from  Table  IX  is  that 
while  all  of  the  optimization  techniques  produced  signifi- 
cant improvements  in  the  accuracy  of  the  maximum  derivative 
guantity  Tmax,  the  same  cannot  be  said  for  the  maximum 
solution  quantity  y^max  and  the  integral  quantity  T.  One 
might  even  conclude  that  grid  optimization  is  not  cost 
effective  in  the  computation  of  these  values  since  the 
uniform  grid  provides  estimates  which  are  nearly  as  accu- 
rate, and  in  some  cases  better,  than  the  optimum  grid.  This 
conclusion  would  be  correct  if  the  solution  maximum  and  its 
integral  were  the  only  resultants  of  interest  in  performing 
the  analysis.  Since  the  purpose  of  this  study  was  to  find  an 
optimum  grid  which  produced  acceptable  errors  for  all  the 
resultants,  the  uniform  grid  is  clearly  inadequate. 
Moreover,  since  in  the  majority  of  engineering  problems  it 
is  the  derivative  of  the  solution  variable  which  is  of 
primary  interest,  it  deserves  special  consideration  in 
making   this   comparison. 

Furthermore,  one  might  conclude  that  the  reason  the 
error  in  the  maximum  solution  variable  is  larger  for  the 
optimal  grid  is  because  the  strain  enargy  density  variation 
criterion  always  concentrates  the  degrees-of -freedom  in  the 
vicinity  of  the  maximum  derivative  value,  which  in  this  case 
does  net  coincide  with  that  of  the  maximum  solution  variable 
value.  However,  such  a  conclusion  is  incorrect  and  might 
erroneously  lead  one  to  attempt  refinement  where  the  maximum 
solution  variable  occurs  in  an  effort  to  improve  its  esti- 
mate. Other  investigations  have  revsaied  that  in  nearly  all 
cases  the  maximum  accuracies  for  ail  three  solution  resul- 
tants are  obtained  by  refinement  ia  the  regions  of  highest 
strain  energy  density  variation  [Ref.  11].  The  more  likely 
source  of  increasing  error  for  the  optimal  grids  is  the 
element  distortion  which  was  encountered  in  ail  but  two 
techniques,        selective    p-version      refinement    and      subdomain 
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isolation.  Such  distortion  can  be  avoided  but  would  require 
more  sophisticated  refinement  techniques  than  ware  available 
for    this   investigation. 

A  reasonable  choice  for  the  optimum  grid  in  Table  IX 
would  be  one  for  which  all  three  values  in  parentheses  are 
as  large  as  reasonably  possible  taking  into  consideration 
the  number  of  cycles  required  to  provide  such  accuracy. 
Based  on  such  a  criterion,  the  author  is  partial  to  subdo- 
main  isolation  for  the  solution  of  two-dimensional  problems 
using  linear  elements,  and  selective  refinement  for  finite 
element  solutions  using  quadratic  elements.  Clearly,  before 
a  concrete  recommendation  could  be  made  for  a  wide  range  of 
applications,  many  more  problems  would  have  to  be  studied, 
but  these  two  techniques  were  fairly  simple  to  implement  for 
a  standard  finite  element  code  and  the  accelerated  conver- 
gence of  the  solution  resultants  of  interest  was  impressive. 
Conceivably,  even  greater  solution  accuracies  might  be 
obtained  by  using  two  or  more  of  these  techniques  in 
combination . 

Here  again  the  crucial  element  in  selecting  the  proper 
optimization  strategy  is  the  precise  definition  of  the 
purpose  for  which  the  finite  element  analysis  is  to  be 
performed.  The  results  of  Table  IX  tend  to  support  the 
following  recommendations  and   conclusions: 

(1)  Regardless  of  the  optimization  strategy  chosen, 
higher  order  elements  are  indispensable  if  high 
accuracies  for  integral  solution  quantities  are 
desired. 

(2)  If  the  maximum  derivative  of  the  solution  variable 
is  of  greatest  concern,  the  strictly  local  refine- 
ments employing  subdomain  isolation  techniques  car- 
provide  exceptional  accuracy  for  a  minimum  number  of 
refinement   cvcles. 
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(3)  If  the  maximum  solution  variabla  value  occurs  at  a 
point  in  the  domain  removed  from  the  vicinity  of  the 
maximum  derivative  value,  then  its  best  estimate 
will  likely  be  obtained  using  a  reasonably  fine 
uniform  grid  and  selectively  subdividing  elements  in 
the  regions  cf  large  strain  energy  density 
variations. 
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VII.    CONCLUSION    AND    RECDHijgNDATIONS 

The  purpose  of  this  paper  has  been  to  present  an  over- 
view of  some  readily  employed  finite  element  grid  optimiza- 
tion methods  and  to  demonstrate  their  effectiveness  in  the 
application  to  some  simple  problems.  This  work  is  by  no 
means  ail  inclusive  and  the  subject  is  still  in  its  infancy. 
While  there  are  many  competing  approaches  to  the  problem, 
there  is  much  more  research  to  be  done  before  any  one 
becomes  widely  accepted  as  a  standard  analysis  tool.  Because 
of  the  limited  time  and  resources  available,  some  of  the 
more  sophisticated  refinement  criteria  and  techniques  which 
have  been  developed  have  not  been  examined  in  any  detail. 
Instead,  the  approach  has  been  to  examine  those  techniques 
which  can  be  easily  incorporated  in  a  basic  finite  element 
code.  However,  it  is  likely  that  some  recently  developed  and 
rather  elaborate  self-adaptive  grid  optimization  codes  will 
soon    be   available. 

Also,  this  paper  has  not  examined  the  important  classes 
of  problems  in  dynamic  and  nonlinear  analysis.  There  is 
considerable  ongoing  research  in  the  extension  of  these 
techniques  to  such  problems,  but  the  increased  complexity  is 
evident.  For  example,  in  vibration  analysis  there  is  an 
optimum  grid  for  each  unique  eigenvalue,  but  it  is  for  these 
types   of    problems   that  grid    optimization    is   most    promising. 

At  the  beginning  of  this  paper  it  was  stated  that  the 
gos,l  of  grid  optimization  was  to  obtain  maximum  solution 
accuracy  for  a.  given  analysis  cost.  Throughout  this  paper  it 
has  been  shown  that,  prior  to  successfully  embarking  upon 
such  a  strategy,  the  underlying  purpose  of  the  analysis  must 
be  explicitly  defined.  Hopefully,  it  has  been  demonstrated 
that    grid      optimization   is    by      no   means    an      unrealistic   goal 
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and    is   far      more   attractive   than   the      non-adaptive    practices 
widely   used  today. 

The   following     are  recommendations      for    future      research 
topics: 

(1)  Investigation  cf  more  sophisticated  refinement 
criteria  based  on  element  residuals  and  reliable 
error   estimates. 

(2)  Investigation  of  grid  optimization  techniques 
employing    adaptive   application    of   the   p-version. 

(3)  Implementation  cf  a  finite  element  preprocessor  for 
performing   hierarchical   grid   refinement. 

(4)  Implementation  of  a  self-adaptive  finite  element 
code. 

(5)  Application  of  grid  optimization  techniques  to  prob- 
lems  of  dynamic  and  nonlinear    analysis. 
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APPENDIX    A 
FORMULATION   0?    THE    ELASTIC    CABLE    PROBLEM 

A.       PBOBLEM    STATEMENT 

Consider  a  perfectly  elastic  cable  initially  stretched 
between  two  fixed  points  a  distance  2L  apart  and  under 
tension    I.      If      the    cable   bears   a   distributed      load    per   unit 


T- 


►  T  — X 


VCX) 


Figure  A. 1        Tension  Cable    Under   Distributed   Loading, 

length   f  (x)    as    shown    in    Figure   A.1,       the   governing    differen- 
tial   equation   for   the  downward    defection    v(x)    is: 


l|£    ♦    f(x)    =    0 


(Eqn.    A.1) 


subject   to   the    essential   boundary   condition: 


v     (x   =    ±L)     =   0 


(Eqn.    A. 2) 
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If   the   distributed  load    is  a    supportive   load    provided   by 
a   Winkler   foundation    cf    modulus   k   such    that 

f(x)     =   -    k    v(x) 

and    if      a   concentrated      load   P      is  applied      at   the      midspan. 
Equation    A.  1  becomes: 


t 

dx: 


:  4^  -  k  v  =  o 


(Eqn.    A. 3) 


subject    to   the    natural  boundary  condition: 


x=0 


+  =     - 


P/2 


(Eqn.    A. 4) 


B-       PROBLEM    SOLUTION 

The   analytical    sclution    of      the    two-point   boundary    value 
prcblem    is: 


p/2 

v  (x)    = [tanhXL   cosh  X  x   -   smh  X  x] 

TX 


(0    <    x    <    L) 
(Eqn.    A. 5) 


where      X    =y  k/T    . 


The  finite  element  solution  is  obtained  by  the  Galerkin 
formulation  using  the  consistent  rather  than  the  lumped 
approximation   for   the  distributed  loading. 
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APPENDIX    B 
FORMULATION   OF   THE    TAPERED    BAR    PROBLEM 

A.       PROBLEM    STATEMENT 

Consider  2.    tapered  bar    of      length   L    and   constant    modulus 

of   elasticity  E    fixed  at      one    end.      Tha    cross-sectional   area 

A(x)       varies  linearly   from    A      at    tha      fixed   end   to    A.     at   the 
*    '  J  o  t 

tip.         Let   the      bar    be  loaded   axially    by      a   concentrated   tip 

load   P,      and  a    distributed    load   for    which    the   intensity   q  (x) 

varies   linearly    from    g      at    the   fixed    end   to   q      at    the   tip  as 

o  t 


UCX) 


p 


Figure    B.1         Tapered  Bar   with   Applied   Loads, 

shewn   in      Figure   B.  1  .        Tha    governing      differential    equation 
for    the   axial   displacement    u(x)    is: 


*& 


A  ^ 

A  ax 


«■   a    =    0 


(0    <    x    <    L) 


(Eqn.    3.1) 
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subject   to   the    essential   boundary   condition: 


u    (x   =    0)     ■    0 


(Eqn.    B.2) 


and   the    natural    boundary   condition: 


du 
o5T 


x=L      EAt 


(Eqn.    B.3) 


B.       PROBLEM    SOLUTION 

Let         a    =    1     -    A^_  /A 

t        o 


and 


0    ■    1    -    qt/q 


o* 


f 

J  X 


For    F  (x)    =   P  ♦  I      q(x)    dx   the   solution   is: 
J  x 


u(x)    -  - 


F(x) 

0  Vx) 


dx 


aA  E 
o 


%<1-    ^)     -    L(l-   k    -   - 

2  a*1  a  q 


Inn-  25.)   +    (1-  %-)x     -  §2.     > 

L  2a  4L 

(0  £  x  <_  L)' 
(Eqn.    8.4) 


The  finite  element  solution  is  obtained  by  the  Galerkin 
formulation  using  the  consistent  rather  than  the  lumped 
approximation   for  the  distributed   loading. 
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APPENDIX    C 
FORMULATION    OF    THE    TORSION    PROBLEM 

A.       PROBLEM    STATEMENT 

Consider  a  solid  circular  shaft  of  radius  "a"  and  shear 
modulus  G,  with  a  circular  groove,  or  keyvay,  of  radius  b 
along   a   generator  of   the      shaft    with    the   cross-section    shown 


Figure  C-1        Cross-section   of   Shaft   with  Keyvay. 

in   Figure   C.1-       An    applied    torque   I    will   produce    an   angle   of 
twist      per    unit      length   9.  The      aquilibrium  condition     is 

satisfied  if  a    torsional   stress      function  \\>    exists   such   that 
the    shear   stress   comccnents   are: 


ZX  9y 


and  t     =    -G8  ^ 

zy  3x 
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The    governing   differential    equation    for    the   torsional    stress 
function         ty  (x,y)     is: 


ill    +    ill 

3x2  3y2 


2       =       0 


(Eqn.    C.1) 


subject   to   the    Dirichlet   condition,       ^  =   0      on   the    boundary 


B.       PROBLEM    SOLUTION 

The   solution   of    Eguation   C.1    [Bef.    22:    pp.    141-143]   is; 


iMx,y)    =   a(x-b2— - —  )    -   ^(x2+y2)    +  %- 


0         9 

x^+yz 


(Eqn.    C.2) 


The  maximum  shear  stress  oc  urs  at  point  A  (Fig.  C.1)  and  is, 


Tmax   =   G9  (2a  -  b) 


(Eqn.  C.3) 


The    applied   torque    computed    from    the   area    integral: 


T    =    2G9    /     ip     dA 

J  A 


G9_ 

4 


(4a4 -8a2 b^b4)  a   +    (2a3b+7ab3  )-Un  a 


(Eqn    C.4) 


where       a     =   arcos    (b/2a)  . 

The    strain    energy  per   unit    length   of    shaft    is: 


U»    =   1_  TT2    dA   =  %T9 

2GA 


(Eqn.    C.5) 


10  1 


The  variational  formulation  of  the  finite  element 
approximation  is  presented  in  detail  in  Chapter  6  of 
Reference   25. 
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